{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TensorFlow version: 2.0.0\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "import datetime, os\n",
    "#hide tf logs\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'  # or any {'0', '1', '2'},\n",
    "#0 (default) shows all, 1 to filter out INFO logs, 2 to additionally filter out WARNING logs, and 3 to additionally filter out ERROR logs\n",
    "import scipy.optimize\n",
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker\n",
    "import time\n",
    "from pyDOE import lhs         #Latin Hypercube Sampling\n",
    "import pandas as pd\n",
    "import seaborn as sns \n",
    "import codecs, json\n",
    "from scipy.optimize import minpack2\n",
    "\n",
    "# generates same random numbers each time\n",
    "np.random.seed(1234)\n",
    "tf.random.set_seed(1234)\n",
    "\n",
    "print(\"TensorFlow version: {}\".format(tf.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Data Prep*\n",
    "\n",
    "Training and Testing data is prepared from the solution file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_1 = np.linspace(-1,1,256)  # 256 points between -1 and 1 [256x1]\n",
    "x_2 = np.linspace(1,-1,256)  # 256 points between 1 and -1 [256x1]\n",
    "\n",
    "X, Y = np.meshgrid(x_1,x_2) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Data\n",
    "\n",
    "We prepare the test data to compare against the solution produced by the PINN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_u_test = np.hstack((X.flatten(order='F')[:,None], Y.flatten(order='F')[:,None]))\n",
    "\n",
    "# Domain bounds\n",
    "lb = np.array([-1, -1]) #lower bound\n",
    "ub = np.array([1, 1])  #upper bound\n",
    "\n",
    "a_1 = 1 \n",
    "a_2 = 1\n",
    "\n",
    "usol = np.sin(a_1 * np.pi * X) * np.sin(a_2 * np.pi * Y) #solution chosen for convinience  \n",
    "\n",
    "u = usol.flatten('F')[:,None] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trainingdata(N_u,N_f):\n",
    "    \n",
    "    leftedge_x = np.hstack((X[:,0][:,None], Y[:,0][:,None]))\n",
    "    leftedge_u = usol[:,0][:,None]\n",
    "    \n",
    "    rightedge_x = np.hstack((X[:,-1][:,None], Y[:,-1][:,None]))\n",
    "    rightedge_u = usol[:,-1][:,None]\n",
    "    \n",
    "    topedge_x = np.hstack((X[0,:][:,None], Y[0,:][:,None]))\n",
    "    topedge_u = usol[0,:][:,None]\n",
    "    \n",
    "    bottomedge_x = np.hstack((X[-1,:][:,None], Y[-1,:][:,None]))\n",
    "    bottomedge_u = usol[-1,:][:,None]\n",
    "    \n",
    "    all_X_u_train = np.vstack([leftedge_x, rightedge_x, bottomedge_x, topedge_x])\n",
    "    all_u_train = np.vstack([leftedge_u, rightedge_u, bottomedge_u, topedge_u])  \n",
    "     \n",
    "    #choose random N_u points for training\n",
    "    idx = np.random.choice(all_X_u_train.shape[0], N_u, replace=False) \n",
    "    \n",
    "    X_u_train = all_X_u_train[idx[0:N_u], :] #choose indices from  set 'idx' (x,t)\n",
    "    u_train = all_u_train[idx[0:N_u],:]      #choose corresponding u\n",
    "    \n",
    "    '''Collocation Points'''\n",
    "\n",
    "    # Latin Hypercube sampling for collocation points \n",
    "    # N_f sets of tuples(x,t)\n",
    "    X_f = lb + (ub-lb)*lhs(2,N_f) \n",
    "    X_f_train = np.vstack((X_f, X_u_train)) # append training points to collocation points \n",
    "    \n",
    "    return X_f_train, X_u_train, u_train \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PINN \n",
    "\n",
    "$W \\in \\mathcal{R}^{n_{l-1}\\times{n_l}}$ \n",
    "\n",
    "Creating sequential layers using the $\\textit{class}$ tf.Module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Sequentialmodel(tf.Module): \n",
    "    def __init__(self, layers, name=None):\n",
    "\n",
    "        self.W = []  #Weights and biases\n",
    "        self.parameters = 0 #total number of parameters\n",
    "\n",
    "        for i in range(len(layers)-1):\n",
    "\n",
    "            input_dim = layers[i]\n",
    "            output_dim = layers[i+1]\n",
    "\n",
    "            #Xavier standard deviation \n",
    "            std_dv = np.sqrt((2.0/(input_dim + output_dim)))\n",
    "\n",
    "            #weights = normal distribution * Xavier standard deviation + 0\n",
    "            w = tf.random.normal([input_dim, output_dim], dtype = 'float64') * std_dv\n",
    "\n",
    "            w = tf.Variable(w, trainable=True, name = 'w' + str(i+1))\n",
    "\n",
    "            b = tf.Variable(tf.cast(tf.zeros([output_dim]), dtype = 'float64'), trainable = True, name = 'b' + str(i+1))\n",
    "\n",
    "            self.W.append(w)\n",
    "            self.W.append(b)\n",
    "\n",
    "            self.parameters +=  input_dim * output_dim + output_dim\n",
    "            \n",
    "        self.X = np.zeros(self.parameters) #store iterates\n",
    "        self.G = np.zeros(self.parameters) #store gradients\n",
    "        self.store = np.zeros((max_iter,2)) #store computed values for plotting\n",
    "        self.iter_counter = 0 # iteration counter for optimizer\n",
    "    \n",
    "    def evaluate(self,x):\n",
    "        \n",
    "        #preprocessing input \n",
    "        x = (x - lb)/(ub - lb) #feature scaling\n",
    "        \n",
    "        a = x\n",
    "        \n",
    "        for i in range(len(layers)-2):\n",
    "            \n",
    "            z = tf.add(tf.matmul(a, self.W[2*i]), self.W[2*i+1])\n",
    "            a = tf.nn.tanh(z)\n",
    "            \n",
    "        a = tf.add(tf.matmul(a, self.W[-2]), self.W[-1]) # For regression, no activation to last layer\n",
    "        \n",
    "        return a\n",
    "    \n",
    "    def get_weights(self):\n",
    "\n",
    "        parameters_1d = []  # [.... W_i,b_i.....  ] 1d array\n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "            \n",
    "            w_1d = tf.reshape(self.W[2*i],[-1])   #flatten weights \n",
    "            b_1d = tf.reshape(self.W[2*i+1],[-1]) #flatten biases\n",
    "            \n",
    "            parameters_1d = tf.concat([parameters_1d, w_1d], 0) #concat weights \n",
    "            parameters_1d = tf.concat([parameters_1d, b_1d], 0) #concat biases\n",
    "        \n",
    "        return parameters_1d\n",
    "        \n",
    "    def set_weights(self,parameters):\n",
    "                \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            shape_w = tf.shape(self.W[2*i]).numpy() # shape of the weight tensor\n",
    "            size_w = tf.size(self.W[2*i]).numpy() #size of the weight tensor \n",
    "            \n",
    "            shape_b = tf.shape(self.W[2*i+1]).numpy() # shape of the bias tensor\n",
    "            size_b = tf.size(self.W[2*i+1]).numpy() #size of the bias tensor \n",
    "                        \n",
    "            pick_w = parameters[0:size_w] #pick the weights \n",
    "            self.W[2*i].assign(tf.reshape(pick_w,shape_w)) # assign  \n",
    "            parameters = np.delete(parameters,np.arange(size_w),0) #delete \n",
    "            \n",
    "            pick_b = parameters[0:size_b] #pick the biases \n",
    "            self.W[2*i+1].assign(tf.reshape(pick_b,shape_b)) # assign \n",
    "            parameters = np.delete(parameters,np.arange(size_b),0) #delete \n",
    "\n",
    "            \n",
    "    def loss_BC(self,x,y):\n",
    "\n",
    "        loss_u = tf.reduce_mean(tf.square(y-self.evaluate(x)))\n",
    "        return loss_u\n",
    "\n",
    "    def loss_PDE(self, x_to_train_f):\n",
    "    \n",
    "        g = tf.Variable(x_to_train_f, dtype = 'float64', trainable = False)\n",
    "\n",
    "        k = 1    \n",
    "\n",
    "        x_1_f = g[:,0:1]\n",
    "        x_2_f = g[:,1:2]\n",
    "\n",
    "        with tf.GradientTape(persistent=True) as tape:\n",
    "\n",
    "            tape.watch(x_1_f)\n",
    "            tape.watch(x_2_f)\n",
    "\n",
    "            g = tf.stack([x_1_f[:,0], x_2_f[:,0]], axis=1)\n",
    "\n",
    "            u = self.evaluate(g)\n",
    "            u_x_1 = tape.gradient(u,x_1_f)\n",
    "            u_x_2 = tape.gradient(u,x_2_f)\n",
    "\n",
    "        u_xx_1 = tape.gradient(u_x_1,x_1_f)\n",
    "        u_xx_2 = tape.gradient(u_x_2,x_2_f)\n",
    "\n",
    "        del tape\n",
    "\n",
    "        q = -( (a_1*np.pi)**2 + (a_2*np.pi)**2 - k**2 ) * np.sin(a_1*np.pi*x_1_f) * np.sin(a_2*np.pi*x_2_f)\n",
    "\n",
    "        f = u_xx_1 + u_xx_2 + k**2 * u - q #residual\n",
    "\n",
    "        loss_f = tf.reduce_mean(tf.square(f))\n",
    "\n",
    "        return loss_f, f\n",
    "    \n",
    "    def loss(self,x,y,g):\n",
    "\n",
    "        loss_u = self.loss_BC(x,y)\n",
    "        loss_f, f = self.loss_PDE(g)\n",
    "\n",
    "        loss = loss_u + loss_f\n",
    "\n",
    "        return loss, loss_u, loss_f \n",
    "    \n",
    "    def optimizerfunc(self,parameters):\n",
    "        \n",
    "        self.set_weights(parameters)\n",
    "       \n",
    "        with tf.GradientTape() as tape:\n",
    "            tape.watch(self.trainable_variables)\n",
    "            \n",
    "            loss_val, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "            \n",
    "        grads = tape.gradient(loss_val,self.trainable_variables)\n",
    "                \n",
    "        del tape\n",
    "        \n",
    "        grads_1d = [ ] #store 1d grads \n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            grads_w_1d = tf.reshape(grads[2*i],[-1]) #flatten weights \n",
    "            grads_b_1d = tf.reshape(grads[2*i+1],[-1]) #flatten biases\n",
    "\n",
    "            grads_1d = tf.concat([grads_1d, grads_w_1d], 0) #concat grad_weights \n",
    "            grads_1d = tf.concat([grads_1d, grads_b_1d], 0) #concat grad_biases\n",
    "        \n",
    "        return loss_val.numpy(), grads_1d.numpy()\n",
    "    \n",
    "    def optimizer_callback(self,parameters):\n",
    "                \n",
    "        loss_value, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "        \n",
    "        u_pred = self.evaluate(X_u_test)\n",
    "        error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)\n",
    "        \n",
    "        alpha, count, f_newval = self.LbfgsInvHessProduct(parameters)\n",
    "        \n",
    "        # if line search fails, it returns None, but difficult to handle while plotting\n",
    "        # so change it to 0\n",
    "        \n",
    "        if alpha == None:\n",
    "            alpha = 0\n",
    "        \n",
    "        tf.print(loss_value, f_newval, (loss_value-f_newval)/loss_value, alpha, count)\n",
    "        \n",
    "    def LbfgsInvHessProduct(self,parameters):\n",
    "\n",
    "        self.iter_counter += 1  #update iteration counter \n",
    "\n",
    "        x_k = parameters  \n",
    "\n",
    "        self.X = np.vstack((x_k.T,self.X)) #stack latest value on top row\n",
    "\n",
    "        old_fval,g_k = self.optimizerfunc(parameters) #obtain grads and loss value\n",
    "\n",
    "        self.G = np.vstack((g_k.T,self.G)) #stack latest grads on top row\n",
    "\n",
    "        n_corrs = min(self.iter_counter, maxcor) #for iterations < maxcor, we will take all available updates\n",
    "        \n",
    "        sk = self.X = self.X[:n_corrs] #select top 'n_corrs' x values, with latest value on top by construction\n",
    "        yk = self.G = self.G[:n_corrs] #select top 'n_corrs' gradient values, with latest value on top by construction \n",
    "\n",
    "        #linear operator B_k_inv    \n",
    "        hess_inv = scipy.optimize.LbfgsInvHessProduct(sk,yk) #instantiate class\n",
    "\n",
    "        p_k = - hess_inv.matvec(g_k) #p_k = -B_k_inv * g_k\n",
    "\n",
    "        gkpk = np.dot(p_k,g_k) #term 1 in report\n",
    "\n",
    "        norm_p_k_sq = (np.linalg.norm(p_k,ord=2))**2 # norm squared\n",
    "               \n",
    "        #store the values\n",
    "        self.store[self.iter_counter-1] = [gkpk,norm_p_k_sq]\n",
    "        \n",
    "        def line_search_armijo(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1):\n",
    "            \n",
    "            \"\"\"\n",
    "            Minimize over alpha, the function ``f(xk+alpha pk)``\n",
    "            \n",
    "            Parameters\n",
    "            ----------\n",
    "            f : callable\n",
    "                Function to be minimized.\n",
    "            xk : array_like\n",
    "                Current point.\n",
    "            pk : array_like\n",
    "                Search direction.\n",
    "            gfk : array_like\n",
    "                Gradient of `f` at point `xk`.\n",
    "            old_fval : float\n",
    "                Value of `f` at point `xk`.\n",
    "            args : tuple, optional\n",
    "                Optional arguments.\n",
    "            c1 : float, optional\n",
    "                Value to control stopping criterion.\n",
    "            alpha0 : scalar, optional\n",
    "                Value of `alpha` at start of the optimization.\n",
    "                \n",
    "            Returns\n",
    "            -------\n",
    "            alpha\n",
    "            f_count\n",
    "            f_val_at_alpha\n",
    "            Notes\n",
    "            -----\n",
    "            Uses the interpolation algorithm (Armijo backtracking) as suggested by\n",
    "            Wright and Nocedal in 'Numerical Optimization', 1999, pp. 56-57\n",
    "            \"\"\"\n",
    "            \n",
    "            xk = np.atleast_1d(xk)\n",
    "            fc = [0]\n",
    "\n",
    "            def phi(alpha1):\n",
    "                fc[0] += 1\n",
    "                return f(xk + alpha1*pk, *args)\n",
    "\n",
    "            if old_fval is None:\n",
    "                phi0 = phi(0.)\n",
    "            else:\n",
    "                phi0 = old_fval  # compute f(xk) -- done in past loop\n",
    "\n",
    "            derphi0 = np.dot(gfk, pk)\n",
    "            if derphi0 > 0:\n",
    "                print('Warning dephi0 :' + str(derphi0))\n",
    "            alpha, phi1 = scalar_search_armijo(phi, phi0, derphi0, c1=c1,\n",
    "                                               alpha0=alpha0)\n",
    "            return alpha, fc[0], phi1\n",
    "\n",
    "\n",
    "        def line_search_BFGS(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1):\n",
    "            \n",
    "            \"\"\"\n",
    "            Compatibility wrapper for `line_search_armijo`\n",
    "            \"\"\"\n",
    "            \n",
    "            r = line_search_armijo(f, xk, pk, gfk, old_fval, args=args, c1=c1,\n",
    "                                   alpha0=alpha0)\n",
    "            return r[0], r[1], 0, r[2]\n",
    "\n",
    "\n",
    "        def scalar_search_armijo(phi, phi0, derphi0, c1=1e-4, alpha0=1, amin=0):\n",
    "            \n",
    "            \"\"\"\n",
    "            Minimize over alpha, the function ``phi(alpha)``.\n",
    "            Uses the interpolation algorithm (Armijo backtracking) as suggested by\n",
    "            Wright and Nocedal in 'Numerical Optimization', 1999, pp. 56-57\n",
    "            alpha > 0 is assumed to be a descent direction.\n",
    "            \n",
    "            Returns\n",
    "            -------\n",
    "            alpha\n",
    "            phi1\n",
    "            \n",
    "            \"\"\"\n",
    "            \n",
    "            phi_a0 = phi(alpha0)\n",
    "            if (phi_a0 <=  phi0 + c1*alpha0*derphi0):\n",
    "                return alpha0, phi_a0\n",
    "\n",
    "            # Otherwise, compute the minimizer of a quadratic interpolant:\n",
    "\n",
    "            alpha1 = -(derphi0) * alpha0**2 / 2.0 / (phi_a0 - phi0 - derphi0 * alpha0)\n",
    "            phi_a1 = phi(alpha1)\n",
    "\n",
    "            if (phi_a1 <= phi0 + c1*alpha1*derphi0):\n",
    "                return alpha1, phi_a1\n",
    "\n",
    "            # Otherwise, loop with cubic interpolation until we find an alpha which\n",
    "            # satisfies the first Wolfe condition (since we are backtracking, we will\n",
    "            # assume that the value of alpha is not too small and satisfies the second\n",
    "            # condition.\n",
    "\n",
    "            while alpha1 > amin:       # we are assuming alpha>0 is a descent direction\n",
    "                factor = alpha0**2 * alpha1**2 * (alpha1-alpha0)\n",
    "                a = alpha0**2 * (phi_a1 - phi0 - derphi0*alpha1) - \\\n",
    "                    alpha1**2 * (phi_a0 - phi0 - derphi0*alpha0)\n",
    "                a = a / factor\n",
    "                b = -alpha0**3 * (phi_a1 - phi0 - derphi0*alpha1) + \\\n",
    "                    alpha1**3 * (phi_a0 - phi0 - derphi0*alpha0)\n",
    "                b = b / factor\n",
    "\n",
    "                alpha2 = (-b + np.sqrt(abs(b**2 - 3 * a * derphi0))) / (3.0*a)\n",
    "                phi_a2 = phi(alpha2)\n",
    "\n",
    "                if (phi_a2 <= phi0 + c1*alpha2*derphi0):\n",
    "                    return alpha2, phi_a2\n",
    "\n",
    "                if (alpha1 - alpha2) > alpha1 / 2.0 or (1 - alpha2/alpha1) < 0.96:\n",
    "                    alpha2 = alpha1 / 2.0\n",
    "\n",
    "                alpha0 = alpha1\n",
    "                alpha1 = alpha2\n",
    "                phi_a0 = phi_a1\n",
    "                phi_a1 = phi_a2\n",
    "\n",
    "            # Failed to find a suitable step length\n",
    "            return None, phi_a1\n",
    "        \n",
    "        # wrapper for line_search_BFGS\n",
    "        def ls_function(x):\n",
    "            val, _ = self.optimizerfunc(x)\n",
    "            return val\n",
    "        \n",
    "        alpha, count, _, f_newval = line_search_BFGS(ls_function, x_k, p_k, g_k, old_fval, args=(), c1=1e-4, alpha0=1)    \n",
    "        \n",
    "        return alpha, count, f_newval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Loss Function*\n",
    "\n",
    "The loss function consists of two parts:\n",
    "1. **loss_BC**: MSE error of boundary losses\n",
    "2. **loss_PDE**: MSE error of collocation points satisfying the PDE\n",
    "\n",
    "**loss** = loss_BC + loss_PDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_u = 400 #Total number of data points for 'u'\n",
    "N_f = 10000 #Total number of collocation points \n",
    "\n",
    "# Training data\n",
    "X_f_train, X_u_train, u_train = trainingdata(N_u,N_f)\n",
    "\n",
    "layers = np.array([2, 50, 50, 50, 1]) #3 hidden layers\n",
    "\n",
    "maxcor = 200 \n",
    "max_iter = 5000\n",
    "\n",
    "PINN = Sequentialmodel(layers)\n",
    "\n",
    "init_params = PINN.get_weights().numpy()\n",
    "\n",
    "start_time = time.time() \n",
    "\n",
    "# train the model with Scipy L-BFGS optimizer\n",
    "results = scipy.optimize.minimize(fun = PINN.optimizerfunc, \n",
    "                                  x0 = init_params, \n",
    "                                  args=(), \n",
    "                                  method='L-BFGS-B', \n",
    "                                  jac= True,        # If jac is True, fun is assumed to return the gradient along with the objective function\n",
    "                                  callback = PINN.optimizer_callback, \n",
    "                                  options = {'disp': None,\n",
    "                                            'maxcor': maxcor, \n",
    "                                            'ftol': 1 * np.finfo(float).eps,  #The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol\n",
    "                                            'gtol': 5e-10, \n",
    "                                            'maxfun':  max_iter*1.5, \n",
    "                                            'maxiter': max_iter,\n",
    "                                            'iprint': -1,   #print update every 50 iterations\n",
    "                                            'maxls': 50})\n",
    "\n",
    "elapsed = time.time() - start_time                \n",
    "print('Training time: %.2f' % (elapsed))\n",
    "\n",
    "print(results)\n",
    "\n",
    "# np.savetxt('values_stiff.txt', PINN.store)\n",
    "\n",
    "PINN.set_weights(results.x)\n",
    "\n",
    "''' Model Accuracy ''' \n",
    "u_pred = PINN.evaluate(X_u_test)\n",
    "\n",
    "error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)        # Relative L2 Norm of the error (Vector)\n",
    "print('Test Error: %.5f'  % (error_vec))\n",
    "\n",
    "u_pred = np.reshape(u_pred,(256,256),order='F') \n",
    "\n",
    "# #Residual plot\n",
    "# loss_f, f_plot = PINN.loss_PDE(X_u_test)\n",
    "# plt.scatter(X_u_test[:,0:1], X_u_test[:,1:2], c=f_plot, cmap = 'jet')\n",
    "# plt.axis('scaled')\n",
    "# plt.colorbar()\n",
    "# plt.savefig('Stiff_Helmholtz_residual.png', dpi = 500)\n",
    "\n",
    "# #plot gkpk\n",
    "# plt.semilogy(PINN.store[:,0])\n",
    "# plt.yscale('symlog')\n",
    "# plt.savefig('gkpk_stiff.png', dpi = 500)\n",
    "\n",
    "# #plot norm_p_k_sq\n",
    "# plt.semilogy(PINN.store[:,1])\n",
    "# plt.yscale('symlog')\n",
    "# plt.savefig('norm_p_k_sq_stiff.png', dpi = 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = np.loadtxt('prove_stiffness/non_stiff.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "s2 = np.loadtxt('prove_stiffness/stiff.txt', comments = \"#\", delimiter = \" \", unpack = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison of actual loss and predicted loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Non-stiff case prediction error:0.0038351804643316583\n",
      "Stiff case prediction error:4.1792986909953895e-05\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABC8klEQVR4nO3deVzVVf748dfhsoOCyCKLCLkloqKiZi5puVVqaZa2TabWOGU5NpXV9G37TZNTzVQ6lVqWs5TamKaWmpoaalpuuG+oqCgIsik73Ht+f3wuiMqqXC7L+/l43Af3nvu5n8/7IPLmfM6mtNYIIYQQ5XGwdwBCCCHqNkkUQgghKiSJQgghRIUkUQghhKiQJAohhBAVcrR3ALbg6+urw8LC7B2GEELUKzt37rygtfa7urxBJoqwsDB27Nhh7zCEEKJeUUqdKqtcbj0JIYSoUL1IFEopD6XUDqXUcHvHIoQQjY1dEoVS6gulVLJSav9V5cOUUkeUUnFKqZdKvTUd+KZ2oxRCCAH266OYD/wT+HdxgVLKBHwMDAYSgO1KqeVAMHAQcK39MIVonAoLC0lISCAvL8/eoQgbcHV1JSQkBCcnpyodb5dEobWOUUqFXVXcE4jTWp8AUEotBO4BPAEPIALIVUqt1Fpbrj6nUupJ4EmA0NBQG0YvRMOXkJBAkyZNCAsLQyll73BEDdJak5qaSkJCAuHh4VX6TF0a9RQMnCn1OgHopbWeAqCUGg9cKCtJAGit5wJzAaKjo2WlQyFuQF5eniSJBkopRfPmzUlJSanyZ+pSoqiQ1nq+vWMQojGRJNFwVfffti6NejoLtCz1OsRaVmsy8jL4+y9/5+t9XyPLrwshhKEutSi2A22VUuEYCWIc8FBtBpCRl8GMLTO4kHOB3Ym7eW/Ie7V5eSGEqJPsNTx2AbAVaK+USlBKTdRaFwFTgB+BQ8A3WusDtRlXmHcY8VPjGRg2kH9u/ydnMs9U/iEhhM189913KKU4fPhwpcd++OGH5OTkXPe15s+fz5QpU8qN46233gLgjTfewN3dneTk5JL3PT09r/u6VXF13e666y4yMjIAmDlzJh06dODhhx8mPz+fQYMGERUVxaJFixg3bhzHjh274evbJVForR/UWgdqrZ201iFa63nW8pVa63Za69Za67ftEZuHswcfDP2AQnMhf9n0F3uEIISwWrBgAX379mXBggWVHnujiaIi7777Lk899VTJa19fX/7+97/b5FplubpuK1euxNvbG4BPPvmEtWvX8tVXX7F7924AYmNjGTt2LH/4wx949913b/j6damPos7o0qILd7a9kxVHVlBgLrB3OELYV/YZuHi0Zh/ZlbfWs7Ky2Lx5M/PmzWPhwoUl5Wazmeeff57IyEg6d+7MrFmzmDlzJufOnWPgwIEMHDgQuPKv/MWLFzN+/HgAVqxYQa9evejatSuDBg3i/PnzFcZx9OhRXFxc8PX1LSmbMGECixYtIi0t7Zrj//GPfxAZGUlkZCQffvghAPHx8XTo0IEnnniCjh07MmTIEHJzc6/9Vmdnc/fdd9OlSxciIyNZtGhRmXULCwvjwoULTJ48mRMnTnDnnXfyt7/9jUceeYTt27cTFRXF8ePH6devH+vWraOoqKjS73dFJFGU453b3+GzEZ+RlnvtD4IQwvaWLVvGsGHDaNeuHc2bN2fnzp0AzJ07l/j4eGJjY9m7dy8PP/wwzz77LEFBQWzYsIENGzZUeN6+ffuybds2du/ezbhx4yr9i3vLli1069btijJPT08mTJjARx99dEX5zp07+fLLL/n111/Ztm0bn332Wclf+ceOHePpp5/mwIEDeHt78+23315zrdWrVxMUFMSePXvYv38/w4YNq7Bus2fPLnlv+vTpfP755/Tr14/Y2Fhat26Ng4MDbdq0Yc+ePRXWsTJ1qTO7TokMiMTJ5MSFnAu08Gxh73CEsB+PlpUfYwMLFixg6tSpAIwbN44FCxbQvXt31q1bx+TJk3F0NH59+fj4VOu8CQkJjB07lsTERAoKCiqddJaYmIif3zUrb/Pss88SFRXF888/X1K2efNmRo0ahYeHBwCjR49m06ZNjBw5kvDwcKKiogDo3r078fHx15yzU6dO/OlPf2L69OkMHz6cfv36VatuZfH39+fcuXN07979us8hLYoKpOel88iSR9h4cqO9QxGiUUlLS2P9+vVMmjSJsLAw3nvvPb755ptqDVsvPVeg9FIkzzzzDFOmTGHfvn3MmTOn0mVK3NzcyjzG29ubhx56iI8//rhK8bi4uJQ8N5lMFBUVcebMGaKiooiKimL27Nm0a9eOXbt20alTJ1599dWSDvQbkZeXh5ub2w2dQxJFBdr5tOPwhcN8uedLe4ciRKOyePFiHn30UU6dOkV8fDxnzpwhPDycTZs2MXjwYObMmVNy3724n6BJkyZcunSp5BwBAQEcOnQIi8XC0qVLS8ozMzMJDg4G4F//+lelsXTo0IG4uLgy33vuueeuiKVfv35899135OTkkJ2dzdKlSytsFbRs2ZLY2FhiY2OZPHky586dw93dnUceeYQXXniBXbt2lVm36jh69CiRkZHX9dlikigq4OPuw+3ht/PD0R/IK5TF0YSoLQsWLGDUqFFXlN13330sWLCASZMmERoaSufOnenSpQtff/01AE8++STDhg0r6fCdMWMGw4cP59ZbbyUwMLDkPG+88Qb3338/3bt3v6KDujz9+/dn9+7dZbZmfH19GTVqFPn5+QB069aN8ePH07NnT3r16sWkSZPo2rVrleu9b98+evbsSVRUFG+++SavvvpqmXWrqvPnz+Pm5kaLFjd2+1w1xBnI0dHRuqZ2uFt2ZBn3LryXJ7o9wdwRc2vknELUdYcOHaJDhw72DqPOmDp1KiNGjGDQoEH2DqVaPvjgA5o2bcrEiROvea+sf2Ol1E6tdfTVx0qLohL3tL+H33X5HZ/t+ox/x/678g8IIRqcV155xWZzNGzJ29ubxx577IbPI6OequDzEZ/Tr2U/OrfojEVbcFCSX4VoTAICAhg5cqS9w6i2xx9/vEbOI7/xqsDJ5MS4TuMwW8ycz6p4co4QQjQ0kiiqyNPZk5hTMfSe15vk7OTKPyCEEA2EJIpq6Bnck7OXzjJp+SR7hyKEELVGEkU19Antw5QeU1hxdAVrT6y1dzhCNGgmk4moqCgiIyO5//77b6gzefz48SxevBiASZMmcfDgwXKP3bhxI7/88ku1r1G8/lJDJImiml6/7XX83P2Yvna6bG4khA25ubkRGxvL/v37cXZ2Zvbs2Ve8f70L3X3++edERESU+/71JoqGTBJFNXm7efNinxfZnbSbDfEVLz4mhKgZ/fr1Iy4ujo0bN9KvXz9GjhxJREQEZrOZF154gR49etC5c2fmzJkDgNaaKVOm0L59ewYNGnTF3hEDBgygeJ7V6tWr6datG126dOGOO+4gPj6e2bNn88EHHxAVFcWmTZtISUnhvvvuo0ePHvTo0YMtW7YAkJqaypAhQ+jYsSOTJk1q0H84yvDY6zC111TCvMLw9/BHay17C4uGb92Aa8tCH4B2T0FRDmy869r3bxpvPPIuwOYxV743aGOVL11UVMSqVasYNmwYALt27WL//v2Eh4czd+5cvLy82L59O/n5+fTp04chQ4awe/dujhw5wsGDBzl//jwRERFMmDDhivOmpKTwxBNPEBMTQ3h4OGlpafj4+DB58mQ8PT1LFvt76KGHmDZtGn379uX06dMMHTqUQ4cO8eabb9K3b19ee+01fvjhB+bNm1flOtU3kiiug5PJicGtBxOXFseFnAv4eVy7sqQQ4sbk5uaWrLbar18/Jk6cyC+//ELPnj1LVnxds2YNe/fuLel/yMzM5NixY8TExPDggw9iMpkICgri9ttvv+b827Zto3///iXnKm8V2nXr1l3Rp3Hx4kWysrKIiYlhyZIlANx99900a9asxupe10iiuE5erl7M2DKDsxfPsmXCFmlViIatohaAo3vF77v6VqsFUay4j+JqxUt4g3GLadasWQwdOvSKY1auXFnt65XHYrGwbds2XF1da+yc9Y30UdyAHkE92JqwlU+2f2LvUIRolIYOHcqnn35KYWEhYKyUmp2dTf/+/Vm0aBFms5nExMQyNzO65ZZbiImJ4eTJk0D5q9AOGTKEWbNmlbwuTl79+/cvWZBw1apVpKen26SOdYEkihvw/K3P07VFV/68/s+cSD9h73CEaHQmTZpEREQE3bp1IzIykt///vcUFRUxatQo2rZtS0REBL/73e/o3bv3NZ/18/Nj7ty5jB49mi5dujB27FgARowYwdKlS0s6s2fOnMmOHTvo3LkzERERJaOvXn/9dWJiYujYsSNLliwhNDS0Vutem2T12Bu0O3E3fb/sSxufNmyftB1nR+daua4QtiSrxzZ8snpsLeoa2JVZd87iRPoJGS4rhGiQpDO7BkzoOoHugd0pshSRlpuGj1v19vAVQoi6TFoUNaRzQGfcnNz466a/ciz1mL3DEUKIGiOJooYopXAxufDJ9k/4vw3/Z+9whBCixkiiqEGtfVoztPVQvj/6PReyG+biYEKIxkcSRQ17qd9LZBdm87ctf7N3KEIIUSMkUdSwXsG9uCP8DubumiutCiFEgyCJwgbevv1twrzCOJhS/pr3QojKfffddyilOHz4cKXHfvjhhze0Z8X8+fOZMmVKuXG89dZbABw5coQBAwYQFRVFhw4dePLJJwFjxnbppUOWL1/OjBkzAGMBwl69etG1a1c2bdrE//73Pzp06MDAgQPZt28f48ePv+64a4MkChvoFdKLFQ+uwMPZg0v5lyr/gBCiTAsWLKBv374sWLCg0mNvNFFU5N133+Wpp54C4Nlnn2XatGnExsZy6NAhnnnmGeDaRDFy5EheeuklAH766Sc6derE7t276devH/PmzeOzzz5jw4YNdOrUiYSEBE6fPm2T2GuCzKOwkeCmwSRcSmDC8gn8Y8g/aOnV0t4hCXFdzmSeIbcot0bP6eboVun/iaysLDZv3syGDRsYMWIEb775JgBms5np06ezevVqHBwceOKJJ9Bac+7cOQYOHIivry8bNmzA09OTrKwsABYvXsz333/P/PnzWbFiBX/5y18oKCigefPmfPXVVwQEBJQbx9GjR3FxccHX1xeAxMREQkJCSt7v1KkTBQUFvPbaa+Tm5rJ582ZefvllcnNz2bFjB5MmTeLFF18seT1q1Cg2b97MxIkTGTlyJO+99x4jRoxg4cKFvPjiizf6rbUJaVHYiMnBhIvJhe+PfM/4ZeMb9KYmQtjCsmXLGDZsGO3ataN58+bs3LkTgLlz5xIfH09sbCx79+7l4Ycf5tlnnyUoKIgNGzaUuQBgaX379mXbtm3s3r2bcePG8e6771Z4/JYtW+jWrVvJ62nTpnH77bdz55138sEHH5CRkYGzszNvvfUWY8eOJTY2tmTdKICoqKgr3nv99deJjo7mq6++4r333gMgOjqaTZs2Xe+3yuakRWFD3YO6M73vdN78+U3mx87n8a6P2zskIarNXq3hBQsWMHXqVADGjRvHggUL6N69O+vWrWPy5Mk4Ohq/vsrbR6I8CQkJjB07lsTERAoKCkr2oyhPYmIifn6X95x5/PHHGTp0KKtXr2bZsmXMmTOHPXv2VLN2V/L39+fcuXM3dA5bqvMtCqXUvUqpz5RSi5RSQ+wdT3W92v9VIvwieHHdi6Rkp9g7HCHqhbS0NNavX8+kSZMICwvjvffe45tvvqlWy7z0HjF5eXklz5955hmmTJnCvn37mDNnzhXvlcXNze2aY4KCgpgwYQLLli3D0dGR/fv3VzmusuTl5eHm5nZD57AluyQKpdQXSqlkpdT+q8qHKaWOKKXilFIvAWitv9NaPwFMBsaWdb66zNHBkc9GfEZ6bjovrqub9x+FqGsWL17Mo48+yqlTp4iPj+fMmTOEh4ezadMmBg8ezJw5cygqKgLK30ciICCAQ4cOYbFYWLp0aUl5ZmYmwcHBAPzrX/+qNJYOHToQFxdX8nr16tUl+18kJSWRmppKcHDwNdevjqNHjxIZGXldn60N9mpRzAeGlS5QSpmAj4E7gQjgQaVURKlDXrW+X+/c2vJW3h/yPuMix5FflG/vcISo8xYsWMCoUaOuKLvvvvtYsGABkyZNIjQ0lM6dO9OlS5eSzYOefPJJhg0bxsCBAwGYMWMGw4cP59ZbbyUwMLDkPG+88Qb3338/3bt3L+mgrkj//v3ZvXt3SWtmzZo1REZG0qVLF4YOHcp7771HixYtGDhwIAcPHiQqKopFixZVq74bNmzg7rvvrtZnapPd9qNQSoUB32utI62vewNvaK2HWl+/bD10hvWxVmu9roLzPQk8CRAaGtr91KlTNoy++grMBRxMOYijgyPtfNrJvhWiTpP9KK40depURowYwaBBg2r83Pn5+dx2221s3ry5pN+lNtTX/SiCgTOlXidYy54BBgFjlFKTy/uw1nqu1jpaax1duuOpWvJSYNcLcLHmV391NjkT6BnIY989JreghKhnXnnlFZvN0Th9+jQzZsyo1SRRXXUpUZRJaz1Ta91daz1Zaz3bphezFMKxj2HfGzY5fYBnAC2btmT2jtkcTJZZ20LUFwEBAYwcOdIm527bti0DBgywyblrSl1KFGeB0uPwQqxltcc9CPz6wfn1YKNbcv8Y8g9cHF2YuGKizK0QQtQLdSlRbAfaKqXClVLOwDhgea1HETgU8pIg48aGu5XnJp+beLHPi2xL2Ma3B7+1yTWEEKIm2Wt47AJgK9BeKZWglJqotS4CpgA/AoeAb7TWB2o9OB9rP07mXptd4vnezxPqFcqcnXOwaIvNriOEEDXBLolCa/2g1jpQa+2ktQ7RWs+zlq/UWrfTWrfWWr9tj9jw7mx8tUGHdjEXRxf+M+o/vH7b6yRlJdnsOkLUZyaTiaioKCIjI7n//vtvqDN5/PjxLF68GIBJkyZx8GD5fYQbN27kl19+qfY1wsLCuHChYW4tUJduPdUNLt4waAuEP2LTy/Rv1Z/gpsGcyTzD+azzNr2WEPWRm5sbsbGx7N+/H2dnZ2bPvnIsS/GEu+r6/PPPiYiIKPf9600UDZkkirK4B0GRbYbClebn4ceY/41h2o/TbH4tIeqzfv36ERcXx8aNG+nXrx8jR44kIiICs9nMCy+8QI8ePejcuTNz5swBQGvNlClTaN++PYMGDSI5ObnkXAMGDGDHjh2AMcu6W7dudOnShTvuuIP4+Hhmz57NBx98QFRUFJs2bSIlJYX77ruPHj160KNHD7Zs2QJAamoqQ4YMoWPHjkyaNKlBD06puwN37enCVoj/D/RfDg62+xZ5Onsy+KbBzI+dz9iOY7nn5ntsdi0hbsSA+QOuKXug4wM81eMpcgpzuOuru655f3zUeMZHjedCzgXGfDPmivc2jt9Y5WsXFRWxatUqhg0zFnPYtWsX+/fvJzw8nLlz5+Ll5cX27dvJz8+nT58+DBkyhN27d3PkyBEOHjzI+fPniYiIYMKECVecNyUlhSeeeIKYmBjCw8NJS0vDx8eHyZMn4+npyfPPPw/AQw89xLRp0+jbty+nT59m6NChHDp0iDfffJO+ffvy2muv8cMPPzBv3rwq16m+kURRltxEOLcKLh0Hr/Y2vdRHwz7i51M/M2nFJKKDogluGmzT6wlRX+Tm5hIVFQUYLYqJEyfyyy+/0LNnz5IVX9esWcPevXtL+h8yMzM5duwYMTExPPjgg5hMJoKCgrj99tuvOf+2bdvo379/ybnKW4V23bp1V/RpXLx4kaysLGJiYliyZAkAd999N82aNauxutc1kijK4t3R+Jp5wOaJoolLExbet5A+X/Rh4vKJrH5ktU2vJ8T1qKgF4O7kXuH7vu6+1WpBFCvuo7iah4dHyXOtNbNmzWLo0KFXHFN6p7kbZbFY2LZtG66urjV2zvpG+ijK0sSaHLKO18rlegT34Jmez5CWm0ZGbkatXFOIhmDo0KF8+umnJau5Hj16lOzsbPr378+iRYswm80kJiaWuZnRLbfcQkxMDCdPngTKX4V2yJAhzJo1q+R1cfLq379/yYKEq1atIj093SZ1rAukRVEWj1Awudts0l1Z3hn0DgeSD5Cam4q3m3etXVeI+mzSpEnEx8fTrVs3tNb4+fnx3XffMWrUKNavX09ERAShoaH07t37ms/6+fkxd+5cRo8ejcViwd/fn7Vr1zJixAjGjBnDsmXLmDVrFjNnzuTpp5+mc+fOFBUV0b9/f2bPns3rr7/Ogw8+SMeOHbn11lsJDQ21w3egdtht9Vhbio6O1sWjGq7bugHg3Az6L6300JpyPus8O87tIDM/k4c6PVRr1xXiarJ6bMNXndVjpUVRnt7/hfxk0BZQtXOHzt/Dn5m/zSTmVAxdWnSho1/HWrmuEEJURPooyuPUBNBQlF1rl1RK8fFdH+NscmbE1yNk61QhRJ0giaI8Dq7w2+9h//+r1cu28WnD16O/Jj4jnud+fK5Wry1EaQ3xtrQwVPffVhJFeRxdwJwHyZtq/dJ3t7ubCV0n8N99/+XrvV/X+vWFcHV1JTU1VZJFA6S1JjU1tVrDfaWPoiLeneHcD8beFErV6qVn3TkLHzcfwpuFU2QpwtGGM8SFuFpISAgJCQmkpMjtz4bI1dWVkJCQKh8vv30q0iwKTi2AS3HQtG2tXtrNyY23Br7FgeQDnM44zU0+N9Xq9UXj5uTkVDJjWQi59VSR5j2Mr2nb7XJ5V0dXsgqyGPCvAfxnz3/sEoMQQkiiqIhPNATcAcp+Da9ugd1o5tqMp1c+TWxSrN3iEEI0XpIoKuLUBLq9D7632C2EJi5N+Pq+r3FQDkxYNkE6F4UQtU4SRWUcnI3RT3bU0b8jrw94nd1Juxm7eCxF5uvbsEUIIa6HdGZXZtdzkHsO7rLdHtpVMbXXVE6kneDQhUOcyjxFa5/Wdo1HCNF4SKKojGMTyLX/VqUOyoGZd87kZPpJ0vPSOZRyiJt9b0bV8rBdIUTjI7eeKuMeDAWpUAdu9yilCGsWhsnBxOD/DOahbx/Coi32DksI0cBJoqiMWzBoM+Qm2DsSwGhZRPpFMurmUSw8sJBHljwifRZCCJuSW0+V8WhpfM0+BZ5hdg2lmKPJkZl3zsSszXy641POXTrH9w99j6ezp71DE0I0QNKiqIxXJ2j1IDh6VH5sLVJK8cndn/DGbW/w86mfeXrl03IbSghhE9KiqEzTm6HDn8Ct6uui1KbXB7xOB78OBDcJ5siFI4R5h+Hm5GbvsIQQDYi0KCrjYAJzAeQl2zuScj3Q8QE6+nckMy+TWz6/hSWHltg7JCFEAyKJoio2jYb9b9g7igp5u3rT1KUpGfkZ3PfNfTz+3eNyK0oIUSMkUVSFawDkJto7ikpF+Eew5/d7GNdxHPP3zOfWebdyPO24vcMSQtRzkiiqwi0Q8uw/6a4qvN28+fq+r5lxxwx2J+3mhbUvkJabZu+whBD1mHRmV4VbEKRssssGRtdDKcX0vtMZ0W4EuUW5nEw/yd7zewnwCKCDXwd7hyeEqGckUVSFWxAUZUNBBrg0s3c0VRbhH4HWmvPZ57l13q0kZSXxdI+nefv2t3F2dLZ3eEKIekJuPVVFi0HQ/o9gKbB3JNWmlKKFZwu+f+h7OgV04v2t79Pyw5b887d/kl+Ub+/whBD1gCSKqmjeA8IfAZOLvSO5bhF+Efw66VcW3rcQf3d/nln1DF/u/pLMvEzZ40IIUaE6f+tJKeUBfAIUABu11l/VfhQmyD4NDi7g7F37l69BYyPHck/7e4g5FYO/pz9xaXG8s/kdmro05e3b3yawSaC9QxRC1DF2aVEopb5QSiUrpfZfVT5MKXVEKRWnlHrJWjwaWKy1fgIYWevBAuhCYy5FvB1ylA24OrkypM0QugR0obVPawothXwZ+yXhH4UzbfU0jqYetXeIQog6xF63nuYDw0oXKKVMwMfAnUAE8KBSKgIIAc5YDzPXYoyXOXmCo2e9mEtRHUopvF29WTZuGb9O+pWugV358NcP6flZT2KTYimyyKq0Qgg73XrSWscopcKuKu4JxGmtTwAopRYC9wAJGMkiFnv2qbgGGDvdNVA9g3uydeJWdp7byTOrniG/KJ995/ex7Mgy/N39GRg+UIbWCtFI1aU+imAutxzASBC9gJnAP5VSdwMryvuwUupJ4EmA0NDQmo/OtUW9mXR3I7oHdeeXib+QU5jD+azzfL3va46lHQOgV3Avnuj2BBO6TpCd9YRoROpSoiiT1jobeLwKx80F5gJER0fX/DAet0BIiavx09ZV7k7uhDcL58BTB9iTtIclh5cwb/c8Jq2YRGJWIhO7TqSFZwtJGEI0AnUpUZwFWpZ6HWItqxtuGg/Ne4KlCBzq0rfNtpxMTkQHRxMdHM1fbv8Lf//l77Rr3o5zl87x373/5WTGSUa0G0H/Vv3xcK5be3YIIWpGXfqNtx1oq5QKx0gQ44CH7BtSKX59wK2FMemuESWK0hyUAy/0eQGA1JxUjqQeYX7sfD7d8Sn+Hv58MfIL+rTsg7ebt30DFULUKHsNj10AbAXaK6USlFITtdZFwBTgR+AQ8I3W+oA94iuTOQ9St0Nukr0jqROauzfn85Gfk/x8Mp/e/Sk5hTkMXzCcexbdw9mLZ2XWtxANiGqIs3Kjo6P1jh07avakKVth7a1wy7/hpkdr9twNQFJWEutPrqeDbwcs2sKBlAN8e+hbxnYcy8h2I/F0kf28hajrlFI7tdbRV5c3znso18PdOpIq50zFxzVSLTxb8FAn405hgbmA2Ttns+nUJpYfWY6DcmBA2ADuanMXU2+ZimMjvXUnRH0l/2Oryj0IHD3g4mF7R1LnOZuc+WzEZ3w49EMWH1zM+vj1rDy2kryiPG5rdRvert609GqJi2P9XTtLiMakSonCut5SrtbaopRqB9wMrNJaF9o0urpEKWjSDi4esnck9YaHswePRT3GY1GPAZB4KRGzNrPy2EpeWvcSAZ4BfDD0A4a1GVbJmYQQ9lTVzuwYwFUpFQysAR7FWIajcfGKgItHwNx48mNNCmwSSEjTEEKbhnJ/xP1kF2Qz/OvhRH4SyR++/wOX8i/ZO0QhRBmqeutJaa1zlFITgU+01u8qpWJtGFfd1PYpCB4B5hwwedk7mnprWNthDGs7jHOXzjFj8wx+O/sb289t58iFIziZnPjx+I9E+EYwtM1QnExO9g5XiEavyolCKdUbeBiYaC0z2SakOqx5L2NPiqJscJZEcaOCmgQx886ZAOQW5pJVkMW5S+d45adXMGszQZ5BRAdH0zukN6M7jKaNTxsclGyhIkRtq2qi+CPwMrBUa31AKXUTsMFmUdVVDiZI3gyWfIh40d7RNChuTm64Obnh5+FH0p+SWHp4KYsOLOK3s7+x/MhycgpzuPfme0m8lMiKoyuI9I9kbMex+Hn42Tt0IRq8as+jUEo5AJ5a64u2CenG2WQeRbGf74GULTA6GRzkr9vasCdpD+5O7rg4urD44GKmr5tesgS6t6s33QK78cldnxDUJAhPZ09Zf0qI63RD8yiUUl8DkzH2g9gONFVKfaS1fq9mw6wH/PvD2eVw6Qh4ybLbtaFLiy4lz5/r/Rx/vOWP/Bz/Mz+d/InjacdZH7+etNw0sgqy+GrfV6w/uZ5ugd0YEzGGlk1bEuEXgcmh8d0pFaKmVPXWU4TW+qJS6mFgFfASsBNofInCp6fx9cKvkijsxEE5MDB8IAPDB5aUFZoLyS7MRmuNm5Mb/9n7H76M/RKAqBZRLB27FGeTMzvO7aBncE8CPAKk5SFEFVU1UTgppZyAe4F/aq0LlVINb+2PqvDtAcoRLmyF1uPtHY2wcjI54W3y5oNhHwBwKuMUe87v4UT6CXIKc8jMy6TIUsTDSx4mqyALJwcnfN19aevTljERY3ig4wM4m5xJzk6mlVcrXJ1c7VwjIeqOqiaKOUA8sAeIUUq1AupsH4VNmVzBpztckhnadVkr71a08m51RZnFYuHbB75l65mtJGUlkZSVxMELBzl04RAJFxPIzMvkjv/cgYNyIMAjgLY+bRkQNoChbYbSyb8TLo4uOJuc7VQjIeznuhcFVEo5Wld8rXNs2pkNkHUSCtLAK9IYLivqPbPFzPns8yw5tIQT6SeIz4hnd+Ju4jPjebXfq9x7873EZ8QzdfVUAjwDaObajN4hvXmhzws0dWlq7/CFqBHldWZXKVEopbyA14H+1qKfgbe01pk1GmUNsXmiMOdD5n5wDTTWgBINktaa5OxkTA4mXEwurD2xlvmx80nOTiYlJ4UT6SdQKFY8uILwZuGsjlvN2uNraebWDD8PPyJ8I+jSogvdArtJS0TUCze6euwXwH7gAevrR4EvgdE1E149Y3KBk/+F5I1w5257RyNsRClFgGdAyevRHUYzusPlH/mVR1eyIX4D7Zq3w0E5kJGXwbG0Y2TmZ5Kem45ZmwHYNnEbTiYnlh5aSnxmPMFNggnwCMDfw5+2zdvSI6iHdKyLOq2qLYpYrXVUZWV1hc1bFAD73oJ9r8Pwo9C0rW2vJeqd/KJ89ifv52DKQYa2GUqRpYiXf3qZJYeWkFWQVXKcr7svax5Zg5PJib/E/IWEiwl4uXoR2jSUXiG96N2yN22atZHhvaJW3Oitp63AC1rrzdbXfYD3tda9azzSGlAriSJtF6zuDt0/hPZTbXst0aDkFOaQeCmRpOwkMnIz6BTQCbPFzP9t+D/2Je/jYt5Fzl46S6GlkJ7BPfnkrk8AeOqHpyiyFNG2eVtGth9Jj6AeBDYJxNNZNoUSNeNGE0UX4N9A8QJH6cBjWuu9NRplDamVRGExw9JA8O0Nty2z7bVEo5NXmMf2c9vJzM+kR1APzNrMtB+ncSrjFAdSDpS0Sh6IeIDpfadjsVgY878xdPDtQKR/JB39O9Letz1tm7WluXtzubUlquSG+ii01nuALkqpptbXF5VSfwTqZKKoFQ4m8LsVkmNAW0AWqxM1yNXJlX6t+l1RtmjMIsC4rbX2xFqOpx0nzDuMAI8ALuZfJMIvgn3J+/jx+I9ojD8Ap/aayqOdHyUlJ4XnfnwOPw8/Wni0IMAzgEDPQO646Q46+BoTR7XWNHWVEVziWtXa4e6q9Z2eAz6s0Wjqm9Bx4OQFBeng0tze0YhGwsXRheHthl9RFkwwKx9eCcCl/EvsPb+XE+knaOPThsAmgWQXZBPgGUBKdgqHUw6TmpuKWZt5m7cxKRO7Enfx5PdP4uHkga+7Ly08WxDYJJCnop+ic0BnkrKS+O3cb3i5eOHj5kO4dzhh3mHSd9JI3Mg8ijNa65Y1HE+NqJVbTwDmAsjcJ8NkRb1j0RaSspJwMRmTCI+nH2fxwcUkZSWReCmRxKxEkrKS+Osdf6WTfydWx63m1Q2vXnEOZ5Mz/x71b7oHdufIhSNsPrMZd0d3PF088XXzxdfdl1tCbqGZWzM71VJU140Ojy1L41zCozSTM5g84My30Ob3xmsh6gEH5UBQk8t/3ES1iCKqRVSZx5otZlp5t+KutneRmZ9Jak4q8RnxxKXF0bpZawrMBfx4/Edm/Tbrms+ufGglLb1asv7ketadWIeXqxfert74uPng4+rDg50exNnkTEZeBs4Ozvh7+ONoupFfS8IWKvwXUUpdouyEoAA3m0RU32TsgZ3PglsghI6xdzRC1DiTgwlfd6OFUJ6Phn3EjEEzyC7IJjM/kwvZF0jNTaVrYFdyCnM4lHKI2KRYsgqyuJh/EbM24+TgRJ+WfVBK8ebPb7Li6AoUqiSRhDQNYd7IeTiZnFhzfA2pOan4e/jT3L05Pq4++Hn40dqnNSZlks56G6swUWitm9RWIPVW8D3g2ATi/yuJQjRaSincndxxd3LHz8OPNj5trnj/0+GfljzXWpOZl0labhohXiEUWYqYHD2Z3iG9Sc5O5kLOBVJyUnBQDuQV5ZFVkMXnuz7n17O/XnHOMO8wFt+/GIDn1z5PfEb8Fa2Vdr7tmNx9MiYHE6uOraLIUoSXqxdNXZrSzLUZwU2Dadm0pSSZKpA23o1ydDX20U5YAkV5xmshRLmUUni7eePt5g0YfR13tb2Lu9reVe5nYsbHkJyTTFJWEhdyLpCem46DcqClV0uKLEVEB0bj7uROem46Zy+e5UDyAZKzkxnTYQwWbeGtn9/i9MXTV5zz1pBbmXnnTBwdHHlq5VMAeLl40cS5CU1cmtAjqAejOozCpEysOLoCdyd3mrk2I8AjgBaeLfD18G00W/NKoqgJQXfDqa/h6CyIeMHe0QjR4Dg7OhPSNISQpiFlvv/ekIq3xtkyYQvpeelk5GWQmZdJel46TV2bEtQkiEJLId4u3py9dJakrCSyC7LJLsgmNSeV6KBotNZMXD4Ri7Zccc6xHcfy2m2vgYaJyyfi4+ZDc/fm+Lr74u3qTa/gXnQP6k6RuYjdSbsJ8QqhffP29XLdL0kUNaHlaDjcDRKWQvtpIJ1xQtQpIV4hhHiVnWQAfnj4h2vKtNaYtZkicxE7nthBZn4mGbkZJOckk5ydTBufNriYXEjNTaXAXMChC4dIz00nMz8Tjeap6Kfwdfcl8VIiIxaOAIxBBCFNQ2jn044nuj9Bv9B+pOels/bEWmPEmLMnHs4eeDh50LZ5W3zdfbFYLJi1mSYuTezWgpHfaDXB0RUGroLs01CYDiY/e0ckhLhBSikclSOODo50Dexa7nGtac3O3+8seV1kKeJS/qWSfptw73CWj1vOqcxTxKXFcST1CMdSj3E64zTnmp0jNimWP67+4zXn/dugv3FH+B1sPbOVZ1Y/g4NywNvVu2To8av9X6VrYFfi0+P5Kf4nTMpE98DuDAgbgJPJqUa/F5IoaoqrP+QlQ24iuEqiEKKxcnRwvGLuiI+7DyPajyjzWK01HXw70L9Vf7ILsskqzCKrIIucghw6B3TGz8P4XfJqv1fJzM8kJSeFlOwUUnNSyczP5OzFs6yKW8VbMW8BMLn7ZHq37F3jieK6J9zVZbU24e5qxz6Dnc9Aj4+h9cTav74QolGxaAv5Rflcyr9EkS7CzdENb1fv6x7JZYsJd+Jq4b+DuE9h+xRwbwmBQ+wdkRCiAXNQDrg5ueHmZNtpbY1jbFdtcXSB25Ybt6F+HgmJP9o7IiGEuGGSKGqaewgM3gSuAfDzPca+FUIIUY/Vi1tPSql7gbuBpsA8rfUa+0ZUCY9QGPorHJkJylGWIRdC1Gs2/+2llPpCKZWslNp/VfkwpdQRpVScUuqlis6htf5Oa/0EMBkYa8t4a4xbC+j4CuhCSN0O5iJ7RySEENelNv7MnQ8MK12glDIBHwN3AhHAg0qpCKVUJ6XU91c9/Et99FXr5+oHJ0/IPgVr+xiztoUQoh6y+a0nrXWMUirsquKeQJzW+gSAUmohcI/W+h1g+FXHooyxXjOAVVrr+nXTv8UQaNIe9rwM3pEQONjeEQkhRLXY68Z5MHCm1OsEa1l5ngEGAWOUUpPLOkAp9aRSaodSakdKSkrNRXqjnDzh9nXGMuSbRkPSRntHJIQQ1VIveli11jO11t211pO11rPLOWau1jpaax3t51fHZka7B8LANeDcDDYMhgu/2TsiIYSoMnuNejoLlN5GNcRa1nA1bQuDt8LpBWByA3M+mFzsHZUQQlTKXi2K7UBbpVS4UsoZGAcst1MstccjGNpMBsxwYr4xGkoIIeq42hgeuwDYCrRXSiUopSZqrYuAKcCPwCHgG631AVvHUic4eRrLe8S+DOsGwNlrlzcWQoi6RBYFtJfUHRBzj7HibKc3IfIVe0ckhGjkylsUsF50ZjdIzaNh8C/gfxvs/TPsnm7viIQQokySKOzJsxUMWAkho8DBCQoy7R2REEJco16s9dSgmZyh7//g4iHIPgmOHY2kIYQQdYS0KOoCBxO4BsHhj2Dv67IulBCiTpFEUVe4NIOCNDj4DqyKhIvH7B2REEIAkijqDqWg/3fQ6wvIOQtrekHcZ/aOSgghJFHUKUpB68eNjY/cQ2H385B5CCxme0cmhGjEpDO7LmoWBYO3GFupmnMgcx84+4JHiL0jE0I0QtKiqKucPCB0NDS9GQ7MgBVtYNtEKMq1d2RCiEZGWhR1naMHdHoDtBlOfGEMne1Z5gK6QghhE9KiqA+8boZ+/4PQsRA3B7Y/bezDXZAByZvtHZ0QooGTFkV90v1DYwht4mrI2A87noKULTA2DwrSwZwHnmH2jlII0cBIoqhP3FrAwB/h4lEwZxlJAsBSCEsDjecPNbxFHoUQ9iWJor5RCrzaQ27i5bLCi+Db+8oyIYSoIZIo6iu3QBhxFBLXQXY8XNhqlBdkgrOXXUMTQjQs0pldnzVpC20nw9GPL5ftf8t+8QghGiRJFPWdUnDLvy6/PvwP+8UihGiQJFE0BCZHGHnq8utjMs9CCFFzJFE0FJ6hMCIOPNvAnlcgfZ8xGqoyWsPXCg5/aPMQhRD1kySKhqRJa+j6N2NOxarOsHF45Z/R1gUHd//JtrEJIeotSRQNTcvRxl7cAElr4GJcJavPyrwLIUTFZHhsQ+TX25itnbQOijIhYy84ekLTttceqy3G145/rt0YhRD1hiSKhsrkAsF3Q1E2/PZ7SPoJhu2CS0fBpzs4eV4+TmZzCyEqILeeGjpHD2gz2ei3+C4YfhoAmx+48nbUvv8HiWvsFqIQom6TFkVj4N8X+i2BE/PgzBIIvQ8yYiE/HUxusO818O4MgUPsHakQog6SRNFYBN9lPAAKs4xbUJvHQGGmUZax136xCSHqNLn11Bg5eYJPN2Ml2psmXi7PPmu/mIQQdZYkisbMtxf0mgsOLhDxMvwYDXkX7B2VEKKOkUTR2CkHGJcHgUMhLwmOfQJF+faOSghRh0iiEAb//uARBvteh/95wOpoyDwEhz+Acz9W/nmLuZKJfUKI+ko6s4VBKRi0Ec4shaw4SFxr7HGx6znj/QfNRuujPBvvBOdm0HdRrYQrhKg9kijEZR6t4OY/Gs8LLoI5+/J752PAtyc4upf92YIMUPLjJERDJP+zRdmcmwJNYeRJiP8vODWBPX82Wg1NI8C/H7gFXD5em0GZ7BauEMJ26kWiUEp5AD8Db2itv7d3PI2KZxhEvgrZp+HUIsiz7svt2AR6/xtC7jFuW6XvMh5CiAbHpp3ZSqkvlFLJSqn9V5UPU0odUUrFKaVeqsKppgPf2CZKUSUeoXBvAow4AbetAGdv2DTKmOmdceDKY7PPSMe2EA2IrVsU84F/Av8uLlBKmYCPgcFAArBdKbUcMAHvXPX5CUAX4CDgauNYRWUcHKBJuPHw7QsnPoem7cFU6p/mawdAG+tLhd4PAQONFocQot6yaaLQWscopcKuKu4JxGmtTwAopRYC92it3wGu2WlHKTUA8AAigFyl1Eqti9fGFnbj4g0dnr/8esQxOP4F5J2HE19A3Gzj4R4CXf4KLQaDWwu7hSuEuH72mEcRDJwp9TrBWlYmrfWftdZ/BL4GPisvSSilnlRK7VBK7UhJSanJeEVVNGkDUX+FW+YZy5bfuQdMHpCfBlt/B1sfg/S9xtyMTfdByi/2jlgIUUX1ojMbQGs9v5L35wJzAaKjo2WDBXtr1hnGZoG50NhpT2twamrsj5G4xujbMLlDyEhjb+/uHxod5CYPMJXxY1mQCRcPG8uOCCFqlT1aFGeBlqVeh1jLRENkcjI2UAoZboyg8u5oDLnt9AaYc+DUQjjzLVw8Ykz0ixkJB94Bc8GV59k0GtbcAmZZXkSI2maPFsV2oK1SKhwjQYwDHrJDHMJeXH2h0+vGA4y1pRRgzjOSx55XYO/r4BZkHBv2CKT+ZhxrzjF25asJJ7+C2BfgnjPgIHNAhCiPrYfHLgC2Au2VUglKqYla6yJgCvAjcAj4Rmt9oKLziAbO0cX45e/sBbevh1sXQOuJRutDmYwE0mKwcWzm4Zq77vbfQ27ilTPQhRDXsPWopwfLKV8JrLTltUU95eAAYeOMR2mJayFhKRRlwaXjEPsiFOVCq/sheLQxk7yyYbjmAihIuzz6KuReiP8KZBCdEBWqN53ZopELHAyjL0BhOljyIWk9FGZA4iqM6TbWRQ09WoFrkNE3crUfIiDr+OUFDj3CjHJLUa1VQ4j6SBKFqD9cmxsPgPvTjRng+98yRlNZisC5udHyOPBXCB5pzCZ3CwAXf2je00gSAEU5xi5/B942XsuEQCEqJIlC1F8eLaHXZ1eW5ZwBkxscnQWUuqV024rLz2OnQ/jvLr/OSzEWO6xoGXUhGjFJFKJhCRpmPCyFRosjNxHyU8C7y+Vjjn0ChRfBLRhyz8IPHaDTm9DpNfvFLUQdprRueHPToqOj9Y4dO+wdhqiLcpON5GFyNkZUrYw0kkqzbnDnTntHJ4RdKaV2aq2jry6XFoVoXNz8jUex21bChsEQ/ijkJsH5DcaMcRcfoy/j7PfQ4U/GZMGalLDMWM6kY1UWTxbCvqRFIUT2GeNWlCUP1g+GgvSrDlBw8zRo9RD4dKuZzu+1/SFlk7EuVnXlJkL2KfC95cbjEKIUaVEIUR4P64oyFjPc8TPknYP8dCNhpG2HjP1w+AM4/A9w9AQ09F0MnjfByf9Adrx1578OcNNj5W8XW1rKJuPrkVnQrAv49696vEuDjK/Xk2SqK223MVosdIztryXqLEkUQhRzMEGzTkCnUoV/ML5kn4ZT/4OLBwAFusiYqJexB5I3Q9FFYzvYXdOMhQ67vm/sIX7yP4DFSDCO7sbCh56tL59+57PG16r+0s9LvvF6Vkf8V8Zy8ZIoGjVJFEJUhUcoRPzp2vLblhtftTb6MxKWgaOHMQNcm+HgO1CYeeVnAodde578NKMFk5MAzj5GK+fQ+xDxojFPpFj2qZqrU1Uc+9RYX0trmW/SiEmiEKImKAUhI4xHafcmQOElKLpk/XoRHFzBwRkubDWG7rYYBNknjUSz740rP+/gAm2ehJTNxvazyZtqq0YGc47xVZuNFpJolORfXghbcvI0HgReWX7bMuOrucC4haUARy+jn6Mg1diH/PQiY+mSnDOw9RFjccTSVrSH5r3g/Hoj4Yw8abRmck4bK++6NL98rKXIGAaszcYtsOpOLrQUwsVD4B5qLN4oGhVJFELYk8n58vOmbYxHsR7/NL6aC+GODZBzzugbOfeDMUor56yxAVTx6rffldooMnQsdHzZeP/EF8YtrWKd/wLtngXnJlWPUxfBys7GUihDf61+PUW9JolCiLrO5HTlUNhWD1x+ri3GqKzj88A1ADL2Gi2RokvGSCyPVtCs+5WJYu+rxgOMhRFvX2dMPjyzxPi8ydWYS5K88fJnLNaNpIr3BalI4UVY1dVYLt635/XWWtQhkiiEqM+Ug7HtbPRHpQoXXn7aeoLxKJaTBEk/GrsIXjpiJARzHmCBtB1wbqXx2lJqJ8HgEXAp7vLrgkxjsuC5H8C9pdHH4ugBLr5GPClbIesErOkFLYYY62yVbjmJekcm3AkhrqS10adhzjX6JvKTjZbLuZXgHgJN2xkTBos7usHou+i/xHj+y+/g4sErzzk62dhLxDO89upRE/LTYG1f6LcYvCLsHY3NyYQ7IUTVKGXc7ire06N4aXevCKOvwlIAPT4x5li4h4BXJ3B0A7cQQEPrSXDqK6MlUjyK66fbIXO/8bzr++B/m/UXrwPkJRkr/prcjfM42ODXUn4qpO+BFrdX73PZp4xO/JRfGkWiKI+0KIQQtpObYiQNzzYQU2rocOj9EDHdaLmsuWopEuUIN42Hdk8bQ4q3PmoME3awbplrcoGg4RD2KJz73uhbaXkvKGdj0mSLwcZclLzzkHvOeL7zWUhaC3fuMXZKbPsUuPpVHn/WSVh+E3T/CNo/W5PfmTqpvBaFJAohRO3JuwCZB4w5IR6twJxvtEzMedZbXfnG1+a9IGAAFGTAnleM4yx5xldznrFUSl6KMaHxan0WQpM2EL8QDr9fdhy+fcA92JqAnIzO/JB7oHkPSNoA51aAdyc4/7N1F0Xg/ovgVI2RYvWQJAohRMOitbHib36K0YluKTJujXmEGhMas04ZSakw3Wh1pO4wJjyCseshGLfRtNnog+nwvDFhcvU1vycNHZ6HiD+Di3etVM8eJFEIIUR5tAa0kTC0xWjVFGYafRRHPoIz314+tvPbEHw3pMfCnj8bLRIHZ+tXJ7j5T8Yqwxn74PgX1vecjZFfTk3hpgnQ9GZAGf0xJrc6szyKdGYLIUR5lALU5RnrJmdjBrpHKPj3g6I8Y+Ji9mnw7WXMUXELMiYgWgqM0WGWAuPh6GHc0jLnGcOELYXGLTVLoTHyy6+/kSBOf2u9daYu9784uEDv+cYosrPL4ORXxrHK0frVyegvcQ8yzmXyMAYAmNyM+S82IolCCCEq4+gK7Z66sswj1FhipTxNWht9KaVpizVxFELgEKPMnFvqkW/cFnN0M5Z0cfW3jjSz3lYzZxlJJ/ecsWjkqVJzZlxbQNunocM0I1nVIEkUQghRW5TD5ZFbAQOMR3naP208ymIxG4tFekXAiX9DVhw0bQ8ewUbHfA2TRCGEEPWNgwkCBhqPyP+z/eVsfgUhhBD1miQKIYQQFZJEIYQQokKSKIQQQlRIEoUQQogKSaIQQghRIUkUQgghKiSJQgghRIUa5KKASqkU4NR1ftwXuFCD4dQHUufGQercONxInVtpra/ZqKNBJooboZTaUdbqiQ2Z1LlxkDo3Draos9x6EkIIUSFJFEIIISokieJac+0dgB1InRsHqXPjUON1lj4KIYQQFZIWhRBCiApJohBCCFEhSRSlKKWGKaWOKKXilFIv2TueG6GU+kIplayU2l+qzEcptVYpdcz6tZm1XCmlZlrrvVcp1a3UZx6zHn9MKfVYWdeqC5RSLZVSG5RSB5VSB5RSU63lDbnOrkqp35RSe6x1ftNaHq6U+tVat0VKKWdruYv1dZz1/bBS53rZWn5EKTXUTlWqMqWUSSm1Wyn1vfV1g66zUipeKbVPKRWrlNphLau9n22ttTyMfhoTcBy4CXAG9gAR9o7rBurTH+gG7C9V9i7wkvX5S8DfrM/vAlYBCrgF+NVa7gOcsH5tZn3ezN51K6e+gUA36/MmwFEgooHXWQGe1udOwK/WunwDjLOWzwb+YH3+FDDb+nwcsMj6PML68+4ChFv/H5jsXb9K6v4c8DXwvfV1g64zEA/4XlVWaz/b0qK4rCcQp7U+obUuABYC99g5puumtY4B0q4qvgf4l/X5v4B7S5X/Wxu2Ad5KqUBgKLBWa52mtU4H1gLDbB78ddBaJ2qtd1mfXwIOAcE07DprrXWW9aWT9aGB24HF1vKr61z8vVgM3KGUUtbyhVrrfK31SSAO4/9DnaSUCgHuBj63vlY08DqXo9Z+tiVRXBYMnCn1OsFa1pAEaK0Trc+TgADr8/LqXi+/J9bbC10x/sJu0HW23oKJBZIx/uMfBzK01kXWQ0rHX1I36/uZQHPqWZ2BD4EXAYv1dXMafp01sEYptVMp9aS1rNZ+th2vN2pRv2mttVKqwY2NVkp5At8Cf9RaXzT+eDQ0xDprrc1AlFLKG1gK3GzfiGxLKTUcSNZa71RKDbBzOLWpr9b6rFLKH1irlDpc+k1b/2xLi+Kys0DLUq9DrGUNyXlrExTr12RreXl1r1ffE6WUE0aS+EprvcRa3KDrXExrnQFsAHpj3Goo/iOwdPwldbO+7wWkUr/q3AcYqZSKx7g9fDvwEQ27zmitz1q/JmP8QdCTWvzZlkRx2XagrXX0hDNGx9dyO8dU05YDxSMdHgOWlSr/nXW0xC1AprVJ+yMwRCnVzDqiYoi1rM6x3neeBxzSWv+j1FsNuc5+1pYESik3YDBG38wGYIz1sKvrXPy9GAOs10Yv53JgnHWEUDjQFvitVipRTVrrl7XWIVrrMIz/o+u11g/TgOuslPJQSjUpfo7xM7mf2vzZtndvfl16YIwWOIpxn/fP9o7nBuuyAEgECjHuRU7EuDf7E3AMWAf4WI9VwMfWeu8DokudZwJGR18c8Li961VBffti3MfdC8RaH3c18Dp3BnZb67wfeM1afhPGL7044H+Ai7Xc1fo6zvr+TaXO9Wfr9+IIcKe961bF+g/g8qinBltna932WB8Hin831ebPtizhIYQQokJy60kIIUSFJFEIIYSokCQKIYQQFZJEIYQQokKSKIQQQlRIEoUQZVBK/WL9GqaUeqiGz/1KWdcSoq6S4bFCVMC6TMTzWuvh1fiMo7687lBZ72dprT1rIDwhaoW0KIQog1KqeFXWGUA/6z4A06yL8L2nlNpuXev/99bjByilNimllgMHrWXfWRdxO1C8kJtSagbgZj3fV6WvZZ1J+55Sar8y9h4YW+rcG5VSi5VSh5VSX1lnoqOUmqGMPTj2KqXer83vkWg8ZFFAISr2EqVaFNZf+Jla6x5KKRdgi1JqjfXYbkCkNpatBpigtU6zLq+xXSn1rdb6JaXUFK11VBnXGg1EAV0AX+tnYqzvdQU6AueALUAfpdQhYBRws9ZaFy/nIURNkxaFENUzBGMdnViMZcybY6wTBPBbqSQB8KxSag+wDWMxtrZUrC+wQGtt1lqfB34GepQ6d4LW2oKxPEkYxpLZecA8pdRoIOcG6yZEmSRRCFE9CnhGax1lfYRrrYtbFNklBxl9G4OA3lrrLhhrMrnewHXzSz03A8X9ID0xNuQZDqy+gfMLUS5JFEJU7BLG1qrFfgT+YF3SHKVUO+uKnlfzAtK11jlKqZsxtqQsVlj8+atsAsZa+0H8MLazLXdFU+veG15a65XANIxbVkLUOOmjEKJiewGz9RbSfIy9D8KAXdYO5RQub0FZ2mpgsrUf4QjG7adic4G9Sqld2lgiu9hSjP0k9mCshPui1jrJmmjK0gRYppRyxWjpPHddNRSiEjI8VgghRIXk1pMQQogKSaIQQghRIUkUQgghKiSJQgghRIUkUQghhKiQJAohhBAVkkQhhBCiQv8fNsWyGfAxQjUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots()\n",
    "plt.semilogy(s1[:,0], alpha = 0.2, color = 'orange', label = 'Actual (Non-stiff)')\n",
    "plt.semilogy(s1[:,1], '--', color = 'orange', label = 'Predicted ')\n",
    "plt.semilogy(s2[:,0], alpha = 0.2, color = 'green', label = 'Actual (Stiff)')\n",
    "plt.semilogy(s2[:,1], '--', color = 'green', label = 'Predicted')\n",
    "# plt.yscale('symlog')\n",
    "plt.xlabel('iterations')\n",
    "plt.ylabel('Loss')\n",
    "plt.legend()\n",
    "# plt.show()\n",
    "# plt.savefig('Compare_actual_and_predicted_loss.png', dpi = 500)\n",
    "\n",
    "# Non-stiff case Error\n",
    "error_1 = np.linalg.norm((s1[:,0]-s1[:,1]),2)/np.linalg.norm(s1[:,0])\n",
    "print('Non-stiff case prediction error:' + str(error_1))\n",
    "\n",
    "# Stiff case Error\n",
    "error_2 = np.linalg.norm((s2[:,0]-s2[:,1]),2)/np.linalg.norm(s2[:,0])\n",
    "print('Stiff case prediction error:' + str(error_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step size plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002073563428166721\n",
      "0.00017585155024532433\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbkUlEQVR4nO3df5RV9Xnv8fdzfgxEo/kBZBkYI9iIC+THaCbRS2svKFW066JtTdCExCArXHsVVqxXS5dNqvbaaEOjITGxpCC5NhUNtUoMIYlerJKgEXQSEaSiJToRI5BiYgzMnHOe+8fe58yZYQbOMGfmsL/781prFnP2nLPPd8+MH5959nd/t7k7IiISlkyjByAiIvWncBcRCZDCXUQkQAp3EZEAKdxFRAKUa/QAAEaOHOljx45t9DBERBJl8+bNe9x9VG9fOyrCfezYsWzatKnRwxARSRQz+3lfX1NbRkQkQAp3EZEAKdxFRAJ0VPTcRSR8nZ2dtLe3s3///kYPJXGGDx9Oc3Mz+Xy+5tco3EVkSLS3t3PccccxduxYzKzRw0kMd2fv3r20t7czbty4ml+ntoyIDIn9+/czYsQIBXs/mRkjRozo9188CncRGTIK9iNzJN83hbuISIAU7iKSGmbGtddeW3m8ZMkSbrzxxrrse/v27UyfPp2WlhYmTJjAggULAGhra2Pt2rWV561Zs4Zbb70VgN27d3PmmWdy+umn88QTT/Dtb3+bCRMmMGPGjAGPR+EuIqkxbNgwHnjgAfbs2VP3fS9atIhrrrmGtrY2tm3bxsKFC4GDw3327NksXrwYgEcffZTJkyfz7LPPcvbZZ7N8+XK+8Y1vsH79+gGPR+EuIqmRy+VYsGABt99++0Ff27lzJ+eccw5Tpkzh3HPP5ZVXXgHg05/+NIsWLWLatGmcfPLJrF69utd979q1i+bm5srjyZMn09HRwec//3nuu+8+WlpauO+++1i5ciVXX301bW1tXH/99Tz00EO0tLRw0003sWHDBubPn89111038GMd8B5ERPrppu88z9bXfl3XfU4cfTx/8z9OO+zzrrrqKqZMmcL111/fbfvChQu5/PLLufzyy1mxYgWLFi3iwQcfBKLg3rBhAy+88AKzZ8/mkksuOWi/11xzDeeccw7Tpk3jvPPOY968ebz73e/m5ptvZtOmTXz1q18FYOXKlQC0tLQc9LX169ezZMkSWltbB/CdiKhyF5FUOf744/nUpz7F0qVLu23fuHEjH//4xwH45Cc/yYYNGypfu/jii8lkMkycOJFf/vKXve533rx5bNu2jY9+9KM89thjnHXWWRw4cGDwDuQwVLmLyJCrpcIeTJ/97Gc544wzmDdvXk3PHzZsWOVzdwfghhtu4Lvf/S4Q9dUBRo8ezRVXXMEVV1zBpEmT2LJlS30H3g+q3EUkdd773vfysY99jOXLl1e2TZs2jVWrVgHwrW99i7PPPvuQ+7jllltoa2urBPu6devo7OwE4PXXX2fv3r2MGTOG4447jt/85jeDcyCHoHAXkVS69tpru82a+cpXvsLdd9/NlClTuOeee/jyl7/cr/394Ac/YNKkSUydOpXzzz+fL37xi5xwwgnMmDGDrVu3Vk6oDhUr/4nRSK2trT6gm3UUOuChq2D6Yhjxe/UbmIjUzbZt25gwYUKjh5FYvX3/zGyzu/d69jWMyv3NV+G5++GVjY0eiYjIUSGMcPdS9G+p2NhxiIgcJcII93Kol0NeRCTlwgh3L3b/V0Qk5cII93LlXlLlLiICoYS7KncRkW7CCPeSTqiKSG1uueUWTjvtNKZMmUJLSwtPPfUUd9xxB2+//XblORdeeCH79u0DYOnSpUyYMIFPfOITHDhwgJkzZw75nPUjEcbyA6rcRaQGGzdu5OGHH+aZZ55h2LBh7Nmzh46ODubMmcPcuXM55phjALot0fu1r32NRx55hObmZp588kmga7mBo1kglXux+78iIr3YtWsXI0eOrKwVM3LkSFavXs1rr73GjBkzKjfJGDt2LHv27OHKK6/k5Zdf5oILLuC2225j7ty5PP3007S0tPDSSy818lAOS5W7iAy97y2G15+r7z5PmAwX3HrIp5x33nncfPPNjB8/npkzZzJnzhwWLVrEl770JdavX8/IkSO7Pf+uu+5i3bp1la+deeaZLFmyhIcffri+Yx8EgVXumi0jIn175zvfyebNm1m2bBmjRo1izpw5lfXVQ6PKXUSG3mEq7MGUzWaZPn0606dPZ/LkyXzzm99s2FgGUyCVu2bLiMjhbd++nRdffLHyuK2tjZNOOqlhy/IOJlXuIpIab731FgsXLmTfvn3kcjk++MEPsmzZMu69915mzZrF6NGj63Jz6qNBGOGu2TIiUoMPfehD/PjHPz5o+8KFC1m4cGHl8c6dO3v9vNzOSYIw2jKuhcNERKqFEe6lQvyvKncREQgm3FW5iyTB0XDntyQ6ku9bGOFeDnWdUBU5ag0fPpy9e/cq4PvJ3dm7dy/Dhw/v1+t0QlVEhkRzczPt7e3s3r270UNJnOHDh9Pc3Nyv14QR7poKKXLUy+fzjBs3rtHDSI1BCXczuxj4Y+B4YLm7/2Aw3qdCyw+IiHRTc8/dzFaY2RtmtqXH9llmtt3MdpjZYgB3f9DdPwNcCcyp75B7ocpdRKSb/pxQXQnMqt5gZlngTuACYCJwmZlNrHrKX8dfH1xx5b7vt/sH/a1ERJKg5nB398eBX/XY/BFgh7u/7O4dwCrgIovcBnzP3Z+p33D7GlzUjnli++u8daAw6G8nInK0G+hUyDHAq1WP2+NtC4GZwCVmdmVvLzSzBWa2ycw2DfjseVy5GyV+16HWjIjIoJxQdfelwNLDPGcZsAygtbV1YBNf4157lhLFkubQiogMtHL/BXBi1ePmeNvQKnWFe2dRM2ZERAYa7k8Dp5jZODNrAi4F1gx8WP0UV+4ZVe4iIkD/pkLeC2wETjWzdjOb7+4F4Grg+8A24H53f35whnoI8fz2LCUKmusuIlJ7z93dL+tj+1pgbd1GdCSqeu4FVe4iIoEsHFbqassUigp3EZEwwr3Sc3dV7iIihBLu5dkyVqKg2TIiIoGEe9VsGVXuIiKhhHvVPHf13EVEAgv3jKZCiogAoYS7q3IXEakWRrhX2jKaLSMiAqGEu6stIyJSLYxwr15+QG0ZEZFAwl1TIUVEugkj3LtNhVRbRkQkjHDXwmEiIt2EEe7lee5afkBEBAgl3FW5i4h0E0a4d7tZh8JdRCSIcI9uCAWG6zZ7IiKEEu5F3SBbRKRaEOFeqpoKqcpdRCSQcPeqVSE7dYWqiEgg4V6Meu5R5a62jIhIGOGuyl1EpJtAwr26cle4i4gEEu5R5Z4zLfkrIgKBhHt5+QGAYqF4iCeKiKRDEOHu1eFeUriLiAQR7tWVeymeOSMikmZBhLu7wl1EpFoQ4d6tcldbRkQkkHD3rhkypYIqdxGRMMJdlbuISDdhhHtVz728QqSISJoFEe7mqtxFRKoFEe6USpTc4k8V7iIiQYS7eZFOcgCU1JYREQkl3Et0xOHuqtxFREIJ9yKdZIGutd1FRNIskHAvVdoy1VerioikVRjhTpFCuedeUuUuIhJGuHuJTstHD3RCVUQknHAvEIV7SW0ZEZEAwt2dLEUKFrVl0GwZEZEQwj1aNKwQt2U0FVJEJIRwj8O8ZNluj0VE0iz54R732IvlE6oKdxGRAMK91D3c1ZYREQkh3OPKvZSJwt0oUSp5I0ckItJwyQ/3cs89E82WyeAUFO4iknLJD/d4tkzJmgDIUqJQKh3qFSIiwUt+uJe6t2UylFS5i0jqJT/c4567x+GepUShqHAXkXRLfrjHlbtnq8JdbRkRSbm6h7uZnWxmy81sdb333aselXtGlbuISG3hbmYrzOwNM9vSY/ssM9tuZjvMbDGAu7/s7vMHY7C9qlTuXSdUi+q5i0jK1Vq5rwRmVW8wsyxwJ3ABMBG4zMwm1nV0tYhny1BVuXcW1ZYRkXSrKdzd/XHgVz02fwTYEVfqHcAq4KI6j+/wylek5lS5i4iUDaTnPgZ4tepxOzDGzEaY2V3A6Wb2V3292MwWmNkmM9u0e/fuIx9Fef328glVK9GpnruIpFyu3jt0973AlTU8bxmwDKC1tfWI07hYLES3xs4OA6K2jCp3EUm7gVTuvwBOrHrcHG8bUoVCdM9Uq5oK2ampkCKScgMJ96eBU8xsnJk1AZcCa+ozrNpVwl09dxGRilqnQt4LbARONbN2M5vv7gXgauD7wDbgfnd/fvCG2rtCZ0c0xngqpGbLiIjU2HN398v62L4WWFvXEfVToRCdULV81cJhOqEqIimX+OUHisWoLZOpnFB1tWVEJPUSH+6dhU4ALKeLmEREyhIf7sX4hGouH1XuOqEqIhJCuJfbMrmuee6dCncRSbnEh3t5KmQmXz0VUm0ZEUm35Id7sdyWicNdyw+IiCQ/3EvFaCpkNpfHMS0/ICJCAOFebsvksnnIZON57mrLiEi6JT7cS5W2TBYsSxbXDbJFJPWCCfd8LgeZrG6zJyJCAOFengqZzebBMvENshXuIpJuAYR7dEI1n6+u3NVzF5F0S3y4lwrltkw+7rmrchcRSX64l7oqd8tkyVmJgi5iEpGUS364l2fL5HJgGXKm2TIiIgGEe1S5N+WjtkzOXLNlRCT1Eh/uXopvs5eJTqhmTVeoiogkPtzLbRkyWbAMeXOt5y4iqZf4cPf4hCqWVeUuIhJLfLiXZ8tElXuWHK5VIUUk9RIf7l7sXrnnTOu5i4gkP9xL1T33qC2jOzGJSNoFEO5FShiYQSZDDqeotoyIpFwg4R4fRly56wpVEUm7sMI9o/XcRUQggHCnVKRkPSp3tWVEJOUSH+7uRZxs9CC+zZ4uYhKRtEt8uB9UuesG2SIiyQ/3XbyP/xw2IXqQyWgqpIgIAYT7qvxs7jrx76MHlcpdbRkRSbfEh3tHoURTttyWyegG2SIiBBDuncUS+WzXVMiMpkKKiIQQ7k5TrvsJVd0gW0TSLvnhXuheuWcpqnIXkdRLfLgfKJbI5yx6kM2T9aJ67iKSeokOd/forkvDKpV7niwFVe4iknqJDvdiyXGnqy1Trtw1FVJEUi7R4d4RnzjNl0+oZnJkvaAlf0Uk9RId7p2FKMSbqit3CnSqcheRlEt0uB9UuWebospdPXcRSbkgwr0pG8+WidsynUXHXQEvIumV6HDvLMSVe1VbJhPfU1XFu4ikWbLDvVy557qmQmYoAq413UUk1RId7gcOqtxz0WOK6ruLSKolOtwrlXvVRUwAOQq6SlVEUi3h4R5Phcx19dwhqtx1IZOIpFnCw71nW6YpeqwlCEQk5RId7h2VnnvXVEiAnFaGFJGUS3a495wtU27LWFFruotIqiU63A95QlWVu4ikWKLDvaOPqZA5tKa7iKRbosO986BVIbtmy+giJhFJs0SHe0ex56qQ0WyZnC5iEpGUS3a4F3r03CtXqBY0z11EUi3R4d7VlilPhay6iEk9dxFJsWSH+0GVezxbxjRbRkTSLdnhXixhBtlML5W7wl1EUizX6AEMxGdnjueqcz6IWRzu3aZCqucuIumV6HDPZIxhmWzXhqrZMqrcRSTNEt2WOUilLaMlf0Uk3cIKd02FFBEBQgv38toypqmQIpJuYYV71c06dIWqiKRZWOFeWRWySKfaMiKSYmGFe7ZryV9V7iKSZkGGe7QqpMJdRNIrrHCvmgpZVFtGRFIssHCPLmjKmSp3EUm3ul+hambHAl8DOoDH3P1b9X6PQ7w5nslrbRkRSb2aKnczW2Fmb5jZlh7bZ5nZdjPbYWaL481/Cqx2988As+s83sOPNRuHu9aWEZEUq7UtsxKYVb3BzLLAncAFwETgMjObCDQDr8ZPK9ZnmP2QzdNkqtxFJN1qCnd3fxz4VY/NHwF2uPvL7t4BrAIuAtqJAv6Q+zezBWa2ycw27d69u/8j70tG4S4iMpATqmPoqtAhCvUxwAPAn5nZ14Hv9PVid1/m7q3u3jpq1KgBDKOHbJ58RssPiEi61f2Eqrv/FphX7/3WLJNnGEUtHCYiqTaQyv0XwIlVj5vjbY2VzZFXW0ZEUm4g4f40cIqZjTOzJuBSYE19hjUA5Z67ZsuISIrVOhXyXmAjcKqZtZvZfHcvAFcD3we2Afe7+/ODN9QaZfOq3EUk9Wrqubv7ZX1sXwusreuIBqoyz13hLiLpFdbyAwCVK1TVlhGR9Aov3LN53UNVRFIvvHDP5KLb7KnnLiIpFl64Z/PktHCYiKRceOGeKbdl1HMXkfQKL9xVuYuIBBzuqtxFJMXCC/dMnqxukC0iKRdeuGdz5Lyg2+yJSKqFF+6ZqC2jyl1E0iy8cM/myXqBTl2hKiIpFmC4N6nnLiKpF164Z3JkXQuHiUi6hRfucVtGC4eJSJqFF+6ZPBlKlAqFRo9ERKRhwgv3bLxEfUnhLiLpFV64Z/LRv6XOxo5DRKSBwgv3bBMApspdRFIswHAvt2VUuYtIeoUX7nFbJlMq4K7pkCKSTuGFezYKd92NSUTSLLxwjyv3Jl2lKiIpFl64xz33HEU6taa7iKRUeOEeV+55rQwpIikWXrjHUyFzaE13EUmvAMO9qy2j9WWkV+7wo6Ww7TuNHonIoMnVe4dmNgv4MpAF/sndb633exxS3JYZn2ln96/3M+LYYTTlwvt/WKq89QbseBR2b4Ph74IzLodjRx75/jZ+FX74OcjkYO6/wsnT6zZUkd50Fktsf/03tL26jx1vvMW4kccypfldTHj/8QzPZwflPa2ec8HNLAv8B/BHQDvwNHCZu2891OtaW1t906ZN9RnEr3fxu69P5x2/e523fRglDKvPnqVB3sEBMuZ0epa8FSl4hgM09Xs/RvS7fowd4If+Ycayi/fxX8zxv+NVTqj3sEUqOoqlSpt4WC7DgULUVchnjb+/ZAp/cnrzEe3XzDa7e2uvX6tzuP834EZ3Pz9+/FcA7v6FXp67AFgQPzwV2F6HIYwE9tRhP0mRpuNN07GCjjd09Trek9x9VG9fqHdbZgzwatXjduDM3p7o7suAZfV8czPb1Nf/xUKUpuNN07GCjjd0Q3G8akaLiASo3uH+C+DEqsfN8TYRERlC9Q73p4FTzGycmTUBlwJr6vweh1LXNk8CpOl403SsoOMN3aAfb11PqAKY2YXAHURTIVe4+y11fQMRETmsuoe7iIg0nk6oiogEKHHhbmazzGy7me0ws8W9fH2Ymd0Xf/0pMxvbgGHWTQ3H+xdmttXMfmZmj5rZSY0YZ70c7nirnvdnZuZmlujpc7Ucr5l9LP4ZP29m/zLUY6ynGn6fP2Bm683s2fh3+sJGjLMezGyFmb1hZlv6+LqZ2dL4e/EzMzujrgNw98R8EPXxXwJOBpqAnwITezznfwF3xZ9fCtzX6HEP8vHOAI6JP//z0I83ft5xwOPAk0Bro8c9yD/fU4BngffEj9/X6HEP8vEuA/48/nwisLPR4x7A8f4hcAawpY+vXwh8DzDgLOCper5/0ir3jwA73P1ld+8AVgEX9XjORcA3489XA+eaWVJXIDjs8br7end/O374JNH006Sq5ecL8LfAbcD+oRzcIKjleD8D3Onu/wXg7m8M8RjrqZbjdeD4+PN3Aa8N4fjqyt0fB351iKdcBPxfjzwJvNvM3l+v909auPd2BeyYvp7j7gXgTWDEkIyu/mo53mrziSqBpDrs8cZ/up7o7t8dyoENklp+vuOB8Wb2IzN7Ml6YL6lqOd4bgblm1g6sBRYOzdAaor//ffdL3VeFlMYws7lAK/DfGz2WwWJmGeBLwKcbPJShlCNqzUwn+qvscTOb7O77GjmoQXQZsNLd/yFeq+oeM5vk7lq/u5+SVrnXcgVs5TlmliP6027vkIyu/mq64tfMZgI3ALPd/cAQjW0wHO54jwMmAY+Z2U6iPuWaBJ9UreXn2w6scfdOd/9PolVXTxmi8dVbLcc7H7gfwN03AsOJFtkK0aBe0Z+0cK/lCtg1wOXx55cA/8/jsxcJdNjjNbPTgX8kCvYk92PhMMfr7m+6+0h3H+vuY4nOMcx29zqtFz3kavl9fpCoasfMRhK1aV4ewjHWUy3H+wpwLoCZTSAK991DOsqhswb4VDxr5izgTXffVbe9N/qM8hGcgb6QqHp5Cbgh3nYz0X/kEP0yfBvYAfwEOLnRYx7k430E+CXQFn+safSYB/N4ezz3MRI8W6bGn68RtaK2As8BlzZ6zIN8vBOBHxHNpGkDzmv0mAdwrPcCu4BOor/A5gNXAldW/WzvjL8Xz9X7d1lXqIqIBChpbRkREamBwl1EJEAKdxGRACncRUQCpHAXEQmQwl1EJEAKdxGRACnc5ahiZu83s1VmtsnM/sPM1sfbm81sTh32f66Z/fPAR3rQfivji9/jnnq/h0h/KNzlaHMP8G/u3uru44FF8fZzidbGHqipRFc+1lv1+KYSrcEu0jAKdzlqmFmWaB2Vfy9vc/fnzOwPiC7Bv8TM2szs5Hh9kofiCv8nZnZqvI974ztx/cTMfm5mf9zjbaYSXdrOIfbxgJn9HzN73MxeiRdmw8wmxNt+ZmbXmdmOeHu38QGnAyf0fH0Nx3+amT0S/8XyOTP7ipl9+Ei/n5JyjV5/QR/6qP4A1hGtlfOPwO/32D4p/jwPPAr8Xvz4QuDu+POtwBfiz/8A+EmP/bcBow6zjxeB/x1//ifA3URL7z4DnB5v/zrwYB/jawOuq359L8e5Fhhd9Xh4PPbTgHcAPwceaPTPQx/J/dB67nK0uQD4fWA2sM7MPunuDwKnAi/Ez7mYKAT/Nb7JVg54wsyGEwX3TfHztgLvKe/YzPLAu9x9t5l9tI99HEO0TPTt8cvywD7gT4GfuvuzVfuuXoXzVOCF+D1GAP/Q4/XduHvPe4POBJ519+fjsTZV7QMzW+HuV/T2DRPpjcJdjiru7sAGYIOZvQeYYmYbiJZDLcRPm0q0ouDy6tfG67q/6O7l2++dQdyCiU0AttWwj83uXow3TQG2xP+2VT11ElG1Xl6K9013L5jZFKL/CZR6vP5wWoj79GY2GnjL3X8UPz4GeNPMZgCzgL+pOkaRXqnnLkcNMzs/rlgxs/cRtVV+CIyl+700dwHnx3dmwswmx/fJnQp8wMyGm9mxRBX87VWvq/TbD7GPyXQP8SnAz4hu+DI+fm4LMLdqX9Xjq36P6tcfTgddt1j7AtENpMvOIOrjn+ruf6lgl1oo3OVocgmwzcx+CjwMfM6ju/G8AIw0sy1mNg1YQfS7uy0+gfmXccU/FXgAeIroxhBfL1e/seqZMn3to2e4TyKqvO8BWs3sOaJ1uXe6e/mmGZXxAfPoHubl13djZmvjCr3sX4A/NLPtRP9z2Ghmd8Rf+3B8TL/t+1sn0p3Wc5dgmNm/Awvcffsg7Pud7v5W/Pl1RL37v673+/Tx3t8A/ifwt8A6d39iKN5Xkk3hLsEws3bgAz4IN1M2s88R3Rauk+hOQX/hyb5frQRO4S4iEiD13EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRA/x+gcV/wyDSbjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots()\n",
    "\n",
    "sns.distplot(s1[:,3],hist=False,\n",
    "                 kde_kws={\"shade\": False},\n",
    "                 norm_hist=False, label='Non-Stiff')\n",
    "\n",
    "sns.distplot(s2[:,3],hist=False,\n",
    "                 kde_kws={\"shade\": False},\n",
    "                 norm_hist=False, label='Stiff')\n",
    "\n",
    "plt.yscale('symlog')\n",
    "ax.set_ylim([0,5e2])\n",
    "# plt.title('Step length vs iterations')\n",
    "plt.xlabel(r'$ Step \\, length: \\alpha_k $')\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig('prove_stiffness/Step_size_stiff_and_non_stiff.png', dpi = 500)\n",
    "\n",
    "# average step size\n",
    "\n",
    "# Non-stiff problem\n",
    "print(np.mean(s1[:,3]))\n",
    "\n",
    "# Stiff problem\n",
    "print(np.mean(s2[:,3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = np.loadtxt('prove_stiffness/Adam/non_stiff_Adam.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "s2 = np.loadtxt('prove_stiffness/Adam/stiff_Adam.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "s3 = np.loadtxt('prove_stiffness/values_non_stiff.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "s4 = np.loadtxt('prove_stiffness/values_stiff.txt', comments = \"#\", delimiter = \" \", unpack = False)\n",
    "\n",
    "s1 = -s1**2\n",
    "s2 = -s2**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-26.386261790602333\n",
      "-2183651.397699857\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEHCAYAAACa4PC5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABTGElEQVR4nO2deZwcxXn3f88cO7P3rvbSSqvVfSMhCWFuI9sYCBZgbN4APjkMwYntxLxxDMGJr5AXQuIDk5jIHxxMAgoEDMHYYBub05ySEJKMEBI6QELnrlZ7HzNT7x/dPVPdU9XTPdfOrJ7v57Ofnanuqnqmp6eerqeeeh4SQoBhGIZhCk1gvAVgGIZhjg9Y4TAMwzBFgRUOwzAMUxRY4TAMwzBFgRUOwzAMUxRY4TAMwzBFITTeApQqzc3NYsaMGeMtBsMwTNmwfv36I0KIFt1xVjgaZsyYgXXr1o23GAzDMGUDEe1xO35cKBwiqgbwbwBGATwjhLhvnEViGIY57ijbNRwi+ikRHSKiLY7y84loGxHtIKIbzeJPAHhICHEtgIuKLizDMAxTvgoHwD0AzpcLiCgI4F8B/AmARQCuIKJFADoAvGeeFi+ijAzDMIxJ2SocIcRzALodxR8AsEMIsVMIMQrgvwFcDGAvDKUDlPFnZhiGKWcm2uA7FamZDGAomqkAfg7gk0T0YwC/0FUmouuIaB0RrTt8+HBhJWUYhjnOmGgKR8dZAE4H8FEA03QnCSHWCCFWCiFWtrRoPfsYhmGYLJhoCmcf7AqlA8D7UK/rMAzDMEVkorlFvwZgLhHNhKF8LgdwG8x1HQAgImtd581CCCCEQEIAwQAVovmJRWwUOLAZmLIcCEy0Zx8mHwghMBYXSAiBUIAQCqbuk+GxOMbiCQQDhIpgIHlsNJZA18AIIqEgAKA6EkQkFMRYPIE9XYMIBQg10RAqw0FUR4wh8NjQGLbsO4a5bTWoDAdRVRFK/ob39QxhNJZAZTiIUJDQUBlGKBjAWDyBN97rQXNNBERAZTiI1rooEgmBV3Z1YzgWx4ymaoQChNa6CCKhIIbH4vj3Z3fiIwtbEQ0HEAkFMbWhEoEA4a0DvXhi8wFcsKQdoSChuSaC+sowjg2N4acv7MJZc5tRWRFEZTiImc3VAIAnthxA18AoTuyoRzgYwNzWGoSCAWzd34sntxzAqvktCBBhXlstwkHCI6/vw96jQ/jgvGYEAwFMaYiitTaK3uEx3P38LnROqsKM5mqcNL2xIN8nlWsCNiJaC2AVgGYABwF8UwhxNxFdAOAHAIIAfgpgG4DVAAjGPpwxAAkhxJcUbV4H4DoA6OzsPGnPHtc9TEqee/J/sGe4Cv/nY+cjGg5m8ckmALERIBBOKZHRASBcBRAB+9YDdR1AbRvwzG2ASBjnfOim8ZP3OCGeEBiLG9c7EgqAiDAwEsOru7rR2VQFAKiJhNBWF00OWDWREJZ21CMUDBj/A4QfPLUdADCntQYAsHhKHWa11OC+V/bgUO8IKkIBtNZGEAoSzpjTjOqKENY8txMAECBCJBxAVUUQl5/ciaODo7j/lXcBGLcHgfCBmZNw2uwmfP+3b6d9hkuWT0XXwAiee/uIrTwYIFx28jQ010Rwx++2245Nqq7A50+foWzvgiXtGI0l8NTWg7by2mgI15w5E+v2HMUL2+19LWyvw/knTFa2d+XpM3D/q+9iNJawlS+aUofzFqvrXLCkHfMn16Yda6wK48ozZirrfObU6fivl9PHp/MWT8bgaAzPO2Q+ZdYkvLLT6WMFtNRG8JlTp9v6iIQD+PNVc9LO9QIRrRdCrNQdL8sZDhH9FMBHABwSQrRL5ecD+D4MU+FPhBC3EtGlAGYA+BchxC+I6GUAyhACQog1ANYAwMqVK7PSxOF9r2AOgCP9H0JHY1U2TZQmQ0eB7p3A1JPSj733GrDrWeD0Lxvvn/+e8f9DNwHH9gEb7gU6TgZmnAm8/ZvUMZFIb+s4Jp4QGBozvPZrzCfvfT1DWLe7G4un1KEiGMSUhihCwQAOHBvG2lffxSmzJiFAhIXtdaiLhpLKAACaaypw2uxmTKquwM9e3G3r66TpjVgwuRb3mYP9xvd6AAChAOFjS9vx5JYDAID+kRhefKcLAFAXDeHxTfuTbRzsHUbfcAwjMWOWcah3BIAxwxAC2H1kEB2Ng9h/bDhZZ0ZzFbr6R9HVP4qeoZSyAYAPzJiEP77fiwO9Q2nXxhowuwdHbcrm7PktGB6N45Vd3Tg2NIZ9Pam6H1nYig17juLo4CicD9bnnzAZT245gN7hMZtC+eiiNuzuGsD2g/2IJYTt2MXLpuD57UcwMBKztTWlIYqZzTX4w44j6BuO2ZTN5R+Yht/88SCGRuNIJFIyLJlaj/mTa/HQ+r0YGLW397nTpuPVXd3Y1TVgqwMYCuXXfzyAXUcGbOWfP30GfvbibgyN2ZXNn509C3c/vwvDY/bdINefPRvPbDtku16Acc986pTpKBRlqXBg7MG5E8C9VoG0B+ejMLzTXiOix2CY1qYg5b0WNcsKykisDAfTeAx4fwPQvgzY8DNg3vnAtl8BzXOBQ1uB4V6gbQkQqrDX2/GU8X/bE8Z5MhvMr2jva0CHQllNYIQQGBw1zD71lWEQER587T3s6xlCS20EVRVBrJrfiopQAD8xZwAW5yxsw7RJlXjwNeO23XnYGGBOm92EldMbsfZVY6C2nlpjcYF1e+xPsEf6R7GnawBH+kds5bXRELoHRrH21ZRDZ1tdFC21EWzZdwyb9h6znX/l6TNwz4u7MTBqH7S+cNYs/HzDXozFE/j5htRP6sITp2B2SzV+8NR29AyO4Z1D/cljFy+bincO9+Oxje/DaVw5fU4z3j82jLGYwA6pzqdP7cSkqgq8srMbsbi90orORvQMjuKVXcaxZ7elvEuXdjRgcDSOl97pSutrweRaPLnlQFp7J0ytx0gsju0H+5FwVJrVUoPX3+1BLJHAjkN9yfKLTpyKnqFRAEir015ficpwELGEwJv7e5PlH1rQiripTBIJYWuvqSaCmmgI8bjApn2p76IiFEBDVRgAkjNVi7poyOzf/jmrKkKoCAXSPn9lRRDhYCBN3lNmNRV0OaAsDec+9+C8BqABwElEVAFgOoDHCi2j88mkLHj3JWDH74C3HgcGjgCv/xcw2A28+wowZj0JOT7XUE/qtVPZOHn5rtTr0cF8SDzuCCHw8s4ubDvQhzfe60EiIdDVP4Lv//Zt/OCp7Vjz3E78xx92Y9eRAfSPxJJPlIf7RrCnaxB7jw7iP17YlWxv2bQGAEDP0KhtZvCJFVMRDQcxOBrD7946lCz/4qrZqKow1gbksePLH56DhqowRmMJbJEGrb86Zy7qomHEEwJC+i5Xn9iO9vooACQHQotQ0BiAXn/3aNrnJ0LaYNY5qQpEhGCA8M7h/rQ6ATLa29OVfg8ECBAwrqFFS00EZNZRLQEEzAHSOXgChh0dsN+1S6bWg4gQIFLXSfaVLnMgYJRv2dcrnW+YAp39yMeFELaH0AAZ5VadzfvsSj5ABAFgUJr9zG6pSdZxyq2S2daWolz13YUKvPZcrjMcFao9OKcIIWJE9GcwZkXfA/BrIcQfVQ041nCyEiJIQFwA8WKvjQ12A1WT9McTceDAJmP20r0TGD4GTF4CJGJAuNI4Z/cLxv/D29Lrx8dS/0MRYMvDwOG3gWhddvL+4YfZ1RsnhBA43D+CnsEx9A6NYUVnI470jyRNUhYDozGbrXxhex227u9F33AMT/5xd7L8Sx+egzt/vwMjsQRi0gD/oQWt2HawD2PxhG0wmN5UjXCQEE8ARwdGk+XRcBDBQPrAGQoGECBC7/AY+oZTgxYZo2Pa+XXRcPK1bjBzKiLAGGidpdYTMkE9AFpDWv/IWPpBs45cjYhAZkNDY+mBQkiql3ZMoaiS8ikGXAAIKj5vY3XY7Cv98xrrYZYMagWWEMJ2zFB4Kbl3H7Er3+S1k5qzHkasOs7zAY3SJRdFCPvDcThY2DlISSocInoKwGTFoZuFEP/rtz0hxCMAHvFwXs5rONZd7JyqF5Sud4BNDwILLwQmn5B+fLAb2Pm0oSD2vwH0mnb4t39t/PezYL/1F4ZX2WFzkXG4V39u30H9sRInFk/gQO8wOhqrsOvIAB593W6FbW9ImbsAYNX8Fjyz7XDaD/ZDC1qwdX8vYokERsbkJ1zzqVw5iBsDimrwEEKkDR5E6YOgUQ7EFZZdgjHQpA1a0hO383wozrfqOOW0BlIyZyuqOjqsQVAn2wFp1pc6Zs0u9H3JRyxnB/JRh6RjhkK0Kw/VTMpWV6QrbDclAcW1CwYoqbTSHwrM/jXDjnr2Z8x8nt2eMkM2VleknZdPSlLhCCHOyaKaag9OwddqnLjeRIViwFwk7D8AQKFwXvn31Ove/enH/dC7F+j3qEjeuD+3vopMIiHwQ4eH06dO6bQpm9NmN+Gld7pwbND+dD6ntQbPbDuMaMjumZg0tThuB2tQVllek4Na2qBLSIj0QddtJqEeUAlCM1sxZE0da6+PSgOw+p5WySn/1/ejPuY6U1G2p5ZDd6yu0pytaGY4qutw8syU9UBApM1I4DLgW4rXcr6Qy/VyU9pML0CpQJC6a679vjUXTkCgZ1CaMYcKO8MpyzUcDck9OOZazeUowlqNE+smUpkfitC7/W1sFNj+VOZqh7fZ12LciMek9ZwMjKU/jZYCQhiL0r/54wGs39ONh9bvRSIh8D/r30s71+neOtd8Onb+sHU/eO2swW2AMM02fQ6TU8B6VE5rS10One1e169iADxxWkNyNqZTBG53um4AdjumU2yq2ZTcXtfAiPaYfUZi/dfPDI06KawHCd3g7TbH0So2LzMzYT8/tVakqaNVuvqZdDEpyRlOJuQ9OES0F6k9OF8C8GuYe3B0azXFoCj65ugeYOP9hrsxAAx2GQv3NW0ABYD3XjX2vWRiy88LK+c4csfvtiOeEFjYXot4Ajh3cRvWPLczTZG82z2I93tSCvKKD3Ri7avvpnkDhQKpzYUyuqdst8HafTFXpLmyEjQzIk37mvHHWPhWxExXDbSZPJYMXeeiIFQzKbPJtyXPrHAw1Y/TZJXqy32wHxxVBYI3jh08lq6MjL70P9QhRXu6a5pxtqLpw11J2K9DkAgJMk1qSlOpXoG6KTx5tqableaLslQ4QogrNOW/AvCrIoujpKAbave/AdS2G8oGSC32d71j/B1njMUTuPP3OwAY+1f6R2L43GnTce9Le5LnbN1vDG4nz2y0KYvVS9vx+Kb9aT9UawAcc6zFkWkTsPatWFiKxbmnIqWIvD/5AtZirvN8vblJO9B4XOi3ygH7LMLpTaXsQ/kJrPbUJh4AtjWtE80FceuzDI4oBnuXxW/AfTb11gHJq8xje/IeodSsCBhTbHlImdJTZRWmecrte9YrCeu7SJXVV4WT3nt+FEtAd38g/b4sdISUslQ45UBeZzhCACO9QLTeeP9WSejUceGZbYfw+rs9aKuLYngsjo8tbbcNDP3mpry3D9rdcS9ZPhWPvL4vzZnDWiTVzWRiCfVMxvn9WgPS6+/2KMu1P3gXm/vm94/Zyi2X4fTz9QNNXHO+m2lKPlQbDbuuuehmaW5yqYY06/pbMyZ59iPXG1F4qVm4ru9IZVbIG4Ja45BCwobKiuSxY0Pp3nUq89gsM/yMs9xeT1mcLJf3MMn97OlKbfxMrklp+tHPcNikVvYUxGngwCZDySy4AOg7kL92SxghBO78/Q7EEgK10RBmt9bg9NlNyQH9YK9h/nLulLZCfrTU2r1twubTpqxYZrVUJ/cdyE/bk+ujyZmMdnHW5/er+8HvO5q+HqZdHjedBlTn6/agaGdQrmacFK21kaTbttPEl+pDM5jCfQZhP5eSxwTs34dcT3bxdranc0l2HrPNPDx60UUrAtpjhvwGchSCbnMxXrtob9bTXVcAGI2rH3jk6zB9UlVSNp0ZUOulpharYLDCKRB5VTg95mL2BJ7ZHOkfwb6jQxgYjeHUmU02b7G+4Ri2H+xDv/Qj+/KH5+BHv9+RFh/KsgjIv9OaSCi5t0IesDsaK5MDXb80UNRFw9o1mcxrNenOBPqBjdCjfFp2WZPRaC69iUUzg9Kcr/wMinPlOrrZfCYvMBkrGKVuncZoT2Oic5FQNyO1jrmZ4WRSsyKdbMb/odHUjVdhusjrFK9VzxnZwd5eXFkuY8VsJCL88X3FNgWdtyDPcMof6/vL6xc5AWOOxRIJw+STEFi35yhejqXWW2Y0VdvOnd5UhZFYwhbuxBrgnU+Hyf0t0hfQWhdRDqZN1ZHkTOaVXSmX1aUd9RnbkcsjYXnjX0oWa6OeznQDZLcwrTzfxxOsZbbSHLEt9FsKU4+7zd/rDKLGDM2iU5KuEijaq6wI2o75MSm5KjDtd5B+vyyaUpeqo3Ws0H7bLufbaapx3zvj1oPO1FcoWOHkGeu+ytsMJ5EwogFMAGKJBF7bfTT5flJ1BWLSVOSchW14autBm5nslJmTcLh/RLvRTaYmEkoqCtljacHkOqW5aFJNhXLGEgykBln5/FnNNdLAkipvr48q2+8zZ03OWcBsa+Oh1n6vfsp382pTN+TmeqxvZ1DhHedkydR617aMenq3bCeNVRXJY26zAXU/BsLxnVgyAOrIH25eXXrchTjUl/KGm1SdWvfxuqY7TTKPKbtRlCf7camjnhkXf4YzkfbhlASzWoyn87w4DQwfAzb9tzrUTJkQTwgcHRxFPCFsygYwZie9kpms2Vx3kfcwzW2rNWNe2dtVPenNaatJzlh2S9F0m2oqkgOCbUFcUlD2tlMzpaeluGUttdJMSRqq5rbWJs+XZSfpv3z+oKWINIOX1jkASAvEabXjxy1ad771GZzrJ6qBrCqSMjFp13BcFGF636k6bh5dMrXJWZFpzpQ+7aJ2aXaB9MX3VF/q66wjkyJ4rzvlYlwrhQtyutHr2qsMZzDdKcpaaiLmMc2syGd5IeEZTp5pqo4gFBjAcK6PDqMDwEv/lh+hxoFjQ2O26LhOd8vmmgrbugkAaZ0lde3qKg2loNvPIV/mBZNTA/+70g+/saoC3Wb8MdlpQGcukuOKyVhPzYC932mNVcp2pjZWJuWUB5wlHanZgbJ9l5mMcxHZKtd7u+mum/eBVqXc2+qi2mM6LNOPaqBLtZNhz4pE0mxmvpc/knUsE3JATcsV3u0j+VEEUck5wXmvp+rZa6YUr3eTmhW81OnRbD38unmpFRtWOAUgGzu0jb3rDIVThiSEQIDIpmwA+5N/i5khUe8BliqLhIIIkNpEGSCymUoIlPajM85L/bjW7bY7GagGv4pQeth2IBU12fl5AoGU8pKrNSVNHWQzEzZI5iOZpIKC35mJd/dj1/bdBlpHH1XSgO6+yJ8iuaaVoR+ni3rymOP9nBa7aVKWo14KX6Pvy+60YcVYc5vj6BSBaqYckpwGvJI0qen6d6mb/gAj30/+ZnKF4rhROET0cQAfA1AH4G4hxG8K1hfUO4E9s/23+RKlaIzE4ni3exBd/aOY21aTLG+sMlLkyhOUxuowegbHXOKLpXt66byMZCpCAY2JLKVWRp0bORW/unAwoDSBRELBZLmscCzTibEmkCq3ZgGAfTCcnJwd2Nuf0pCaEckPLMl0v47y5BO5zqTmtobjd7Hc8dnsay4ak5rjffLzaXsxjqkcKYD062UpaAt5/SRlmnL7TLBdiJbaiLIfL2RSol6P1UrOE77bcnxW3QzHaYosJmW7hkNEPyWiQ0S0xVF+PhFtI6IdRHSjVS6EeFQIcS2A6wFcVmDZihu8c5x5t3sQG97twZH+UQgAB3tTP/y2umjajV1dETJ/OPZrZCkLub5Rrt+pLzOpusJ1IR4AehVuyOrz1eVuG7GddazI0eSYoaWMR/YKdZLikte25MFdVrwnz5iUbFA2DbVJCk1WBvITfDYzHN1nc9v4KZOadfifQQDp10tuj8jusRgKuu+bsY7JolvekZkUohfZvNVKP5L87rJYd3F+1uR+I5DtAWl2i7vTSiEpW4UDI7/N+XKBlPXzTwAsAnAFES1y1PuGeU5B8e00IARwdHfJBrx0QzYXBQi2ndjW05RMNBxUmnUshSMvvALWWkR6v6rBSf9DVZPJ88lLnxZOZZTKC0PoHZJz0qj7bpK8jY5IT+vTmy3PJXKYtUJJWa2NsIBsmrPPDKskV2G/Jl/tZyO1ec6SS8bqMruhOf16WWm4jXq6NQ/3vuTLIO9n8SqDp358HLN+A9n04yS5Ful44LGyhrJJzQdCiOeIaIajOJn1EwCIyMr6+SYZd9GtAJ4QQmxQtZmPBGxGOz5/0MO9wEumDmyem3W/xUQIgd1dg2k/AOdNHCS9Cki7RC6Dv/yk7txjYTvX5w9VP1D5U1ypOunf+2gs4WnGG0gqKE37sF+H2a2WycReY0aT2rX2BMuV2dGu7Nygw23Tpdc1HJ2iVZ1jMbO5Wn2iqp7vPTWUdCYBgOpISoF7lc+TbI73HZIp0ClfRQ7pAdzuWfnBo7kme9NhrpTzDEeFKuvnVPP1lwGcA+BSIrpeVVkIsUYIsVIIsbKlpSVrIfz43QNIBd8EgCPb9eeVCAICm98/hgO9w7Y0yPWV4bS7WLtx0KeykAe1qQ2pJ3gPzWrPde1T104WP1JdBk29J5JaFud1SIaDcdS31i8A+/kpM5hdSURCqXLfOMx21qK3Ja9Msn8fZqHpTXJ76feWJIan9rwc87NOkmlGYhyzH1wwuU46lk/Z9P3b9p8lk6wVX+OU7AynAFk/7wBwR86CecBps3el6x0j+nMZ8bIUTiYUoGSsrcaqsC2WlIVO31hX6HD1XCyYXKu9/Z3RbpNhUJRKy58GyadiUTksuOFVOSadEhwVIiH1OsUkyTvO5mIeTT3By0qio9HDYr7b7NPhmi4f01bSHlIP6LnIl4UYnislZ5k+1mkyRQZway+7KAj2pH0VmvumGJSswin3rJ+eFU73zoLKkm/kfSBN1RU4KmULrI2GQaRKzqaxd5iMBavQEAm5Kgv5csopglU9KZvQPjH6NKlRuvtzpj50aM/3aP4LSOsotnLrqRv2/T/yQvqAFP5fp9Dc+k6WOyr1yJlQPc7gHB3ZkFN2Z2Xyy0JL+VEec1trPchmfy8rHL+zyuxmOPb3YY27trwmVigmmkmtRLJ+6u3a5c76PaloAc5QK5XhoNrMpZqJALZRe1ZLtee1Czczhm8TmabcDb+OCfp2fJbnqIh07biZKDP3bVdq1hO/dUxXR9uP432zPDi71PPeYuYjfsxjzTXqPVVu7dnNot7qpFzgvcuma0t3fntDVH1iHilbhWNm/XwJwHwi2ktE1wghYgCsrJ9bATw4Hlk/fc1w9q4rqCyFxHkjyzHIMtclWFti9tUtM5SVyw9DDrdinefHhOX3x+jelrrcr0ktoPGxztdTr1ezTNiLC7HHvm2eYz5nlSrqq1JRH4rpPaat43if9BR0k01zvd1kcJYv72zM3I/v+6D4lKxJLROlnPXTWMMZTwnyi4DA4Ggc1RXpt4uXj5nphx0PRDCpugLDihwoQPoPI9NTnsrkpR8w/f/ssnky9tWOv2b07XvsILVfw39bzsFMdljIZqBzfh8VXk1qHvuSN+N6dSSRPeWcx7xcO+dBOcyT/iuyH/Hiyux7Zuwod0ZpLwRlq3BKmYmy8TOeEOgZGk1mz2yUnjZrIyHPayieyzz+MEJBd68q5654t7atY36+rrwttvqcdTn3J6XO9zewa8XPZkDXfDeqY5nKVSJ48URz78t+YJ4UBSObGZjzOgQzmC+NOi7HPN4DyT1M2VwErVz28704M+QKK5wCQMgxlto4ISAQiwuEgwEc6B3GriP2eG5HpQXhusow4or4PZ49x8jDOclTfQ6ohLSpl5+wLZnQL1D7+8H3DaujHuiaGYurZfT/ZOvN1CYHK82nF5guYoBbP5kPepNBnkV7VchWyJssRcjSDGbHy/4cv7N4Z7GVRruQsMIpAOVmUntjb48tf4wXWmsjtj04KbJ//Pf7RKwLM+PXg0yjP9zPV5V7bwKAOo2y0Y7PJ1WfUxavCmqGl02XWayR6BbEjebcHgz8i+Hsy7bp0qN8nfLeIm0/Lt6Lmjru7dnfy4FgtW3leF9Gcth06hVWOAUio0ktNgrExieMzeBoDL3DMfQOjSFA5FvZAABIN5vxXF1RZi+d0hDVnmv0pT7iN4uhX0WRTRgVX+27DIR+gm7qFLLXgcnKWGnU8TYryob6KtlNOLs2vNZrlPvy6FThOOib7GZt9vKGqhz27ni8D3ROLPmEFU4B0EU3tvH6vUD/4aLIE0sk0Dccw0gskWYmk1kwuRZvHehLvnd78A/6GBmsM3uiHaiMHTPKFPWdRckfkO+BM91E5vQgsy0eS4e8fKx8zXD07ednhuNXSTjPl92d/Zrt3PpxsqyjIWM/2fblLK+RY/t5vK9s+2ayMM1mM2tzyuDMJ2XRUJXZHFlKCdjK1i3aL0RUTUTriGh1wfuChzWcAisbAYGXdnbhpZ1deG33Ubx1oM+mbKzcJDIhx03dWK1/qgoHA96fRs3z3mq9AIPhxkynKcrVR1RJ2QBvM6/58o5428J05g+lO0Ne48oFv9fB61pEstyjwrQnq/PXh9tB9yd+t0NZmJOcT/Hyw4VHEUIBr5tPU6/ldSk3cn14mVSdeZ9SNs4bhaIsFY7f1AQmXwfwYHHk0w+GAOyx0/LIq7sMBfPKri5b+BnZu8xCdl+1cD7hZtp5nDYATDvZk6kMAI5N+0jG/hc60gQ70V1i+fRUXnk7qmsCGBlGVchRrwv9Q/U7QPgu99hvU3VqsTybNUmvirNbilaRjVLxJZP8cOHxAoU8uDED9hxJp8yaJPXjS0R32RycKM8OCzzzzgdlqXDgMzUBEX0UwJsADqEIZAzeuev5jG281z2Il3Z22W5iJ0II9A6PYV/PkHGueaqzihws0A/yYq4bL3deh23N5wIzz/Z8d4twZVqZ02pQHcmQ393DQJvcRe84WV4Ilk1Hi9rV18o+6KRey6Y5mRVW0jQHVlIsJ3JKZP+mDn+2e+3Mw3FA/v51eYSyyW3jLJ6sMW+my+dyTKroxcwEQOllafRjrySvbXhVBFZEZqM9PbluyrSZg3NcwykGZbmG4zc1AYBVAKphKKIhIvqVECLtbstnegKn08CmvT3oHYrhzBk1mlopXtrZlXy95f1jtqcYwHiSetWRKlmmtTZiy37oRGUPbquL5DTMHa2aAVBAedPLJdXmwCo89NaaTCTmd02DYK0+NWlCj+jalKMdy8xsSX1vcn76aZPSFSegX+PaeVi9hpbMXYLcZyap833a7h3FXq67VwVh2zzpOM82e/TYl9P8G5PMf9ZDhlFH36J2o7FHs5lnPLZnC1TqsZ+obKnQ1vF5HxSQslQ4GlSpCU4BACHEzQBARFcCOKJSNuZ5awCsAYCVK1dm7dhMMLyJRmJxbN57DAva6/C7rcbkal7dGFo19fYeHcR7R+3BLy0PsoQQ6Oofxb6eIQxJmQ0joQASQtj2aGQKsVKrMJUFiNIHHBD6K1pRM6qZGCq6kYtaalIJxSyMmUWvZp1F/YTtZQCWw9jL+M39MUmzbiXvSZGJafbGyDMWL0j7JXNee8l0INcn3sWy95rLeTpnDOf3LLvjet2P1eb4PmRPy8Vm3h8/8tn7canjckzG9gDh0WnA5hXosafWWvVMyh4gVNM3z3AM8p2awEIIcU/WQvmAyAjZf++Le9A/ErPNRvZ0DWoVjlPZWBzpH8H2Q/3KYys6G20zokXtdTabuMxwqA7RWC+apZvUolHldknAaKgaUDen/klIhaOtS4EVZwHv/yhZZj2Yenl6joSCtjppXUltdGrysKgyjgJ6r5+oYm0LANpq1QpnnuR8IDOnRT2TXdbZoCxPpouG2wCRmrmllyvOVzfju9yJ142a8jG3xW0v0QmcFee1qa874DB/evxQcnvZ7HVxYjORerxGKyUzbK7rPifJbenO999FzpSkwinn1ASAMUAmhEC/GQJe3uBXq/FecfNqk5XN0o56bNp7THtuKKh/NrKyyqiOR8KBNEcHApCAvyf1IfNJ8+3mj+Ijy84C6tUzj4Cv2z3zudXSrE2e4aVC7+f286qKqK9Dk4snn4ppjWoTnH0neUpWeTYhs8jDLCOb1Ate6LCZHb3NSJZObZD60bftdXYxt1Vvmvay18aJPQunG/7vI681dA872fQje9bJKUVs54+DxilXpwEVJZGaADAGPHkh2qJxcDeCPbvSyhNC4OVd+jUZi47Gyow3b5XrTWsolGg4CFQ22o5EQyqvNfe+5AHKGjBjCQEBQnfVTGUk3XQTV/brBHJ5naTIhyTzSr72sskRfmV0g0RQ43AR0ziByCYheYDQOSXYgmR6MA15mZV4vVTNHsw1zmMkmwzzMINwG5yDHmc48iG3TKW54tVMqDNBem27R3LskMcfea3Kfj/xGo4nzNQEqwA0E9FeAN8UQtxNRFZqgiCAn45HagJDPnX5/CO/QWKoAmirRf9IDJv36WcqKtrro9p4Wqm+7dkxK4IBwPIIMw/o4jJlesLdOemDmNX9nPKYyqNNZc4KkCWd95vdy5k6N+dcZzbZonMplxe0ZWSFMDymjvzgRfFWeTDl2GYKHoJZOrHtS/FWJatvW/ed+ukrm28/367ZhXSLlpGz7U7VzKS9rO0UkrJUOKWcmgBQL9pHx3oAAN0Do3jrQG9WmwRDgQDG4vowNAdrFgI4ZDONNVan/2iNkO+qHdHp749Fp6J5cAcAoL+iWXu+6kdqeWrJRyrD1qxH9ynS0f34uvpTi0tVitQJpYjuyVyeYcqfVpcUa2G7eoF5iuyhJZWfOC21kC7fe7YnbI/Ds23zpEe3aK8DqHzaLMc6mNfI0QEvuzsd2LwKszT5eamzwLHmNxyTZ+PePp8X5IdKua3Z0jUdj0exiWRSKx0U3+Sy/ak9p7nsSHe7SXoj7QCAY9LUurGqIvkrts19PESzPrr8z3G4el7yfTzgWK+QF3EVi+fWD1/+UTkRyBwWXf7MkzXeYuWCzjRn31ia+sQ6b6doWO3WJnvryWOopegB/1Gq08/zNjDa1xQ8Khzptd+Edur2vNn8wiGPZrgsZHKrckwaC7JxdrD1I72u1sx0qyNyefFVDiucQjDOkaLTNouGjEH6SPWcVFl1S3pFAl6ZdjW2N30Ym9s+DlRU2e7WkVAd3mxdLZ+epKUm3fNNKw/sPwLZQ0uFfK6bd1I5o3vA9mLXl4vd9ruo2qnWOEN4xesajk7RptdJVXLTUV4DTXodU+UZslclOqNZ7RDjVmuRwwlE/mUENJuLs8E+u0yV29bypPNzvQ+8wgqnABRU33i8D9+ZdDYAcwE1EMKrHVfh3foPpE5YdBHebP0YBsOTMBa0ctoDgkLoqp6DgUir/YasbQAA9EanJMvkzynnsndyrGoG+iJttjLr6TUeqHDdpGrIlZKkVeHSXep4GfR0A4RuA2lQs/s9qDEn6WaGDTlGapa/mxOk/S/OYzo39PT2pNd5meG4oJnle3VqmNvq7eHHruDtpl+doaEYkw83uQpFeRi+mSRuC6L9EWOHz862j6I72I7Z3c8mvZkSAcdaTiiC3uhUbGq/FADwQTyc3pd5R25uuwTtba1Aj/3X0d3vjIGl/vVsavyo9M48JxDA7sbT0RPtwHRpsVOF/MMow7x2vgloFIj8fXsrl0xGGq+5Zo+RkHXI383C9lrtsWzaqwjZG8hq/URq0OmEcKRfvcHMfYaTOtri8eEnGzNhofSNLgW4c22pUPAMpwAUI9vn/toT0nazj4TqgA/dhJ7q2ckyr0+WQPqTnfVuINKCqW3NaecnJAXjd3AhEA7UnoDhcIMt7W8m2urLb4bj93bQ7cqXyyOSk4HsDlshma7kUPmy15w80NhifuU4yjkH0+wW2FO1ZjZ7vy/07Umv8zFj8uZbYEMOheScsWpzN+WorGVsWU4191PnJA/J9vLAcaFwiChARLcQ0Y+I6POF7q82mp07pxesG+ZQzULt05KxId24kfOx8AoAdU4X5/qpdpOJz37k0+WoxJnOrfC4FnA88X5PKkKFPKC83zOsLNeRzZ0yIjmEOMMCuTmL6OgZSs06vAaPdUP+2CGP7Xm9lb3OcOQ9Mc4tCVqTmvRt6KJTONFtmTigzMw7PpTlrzeL9AQXw4g8MAYjxlpBqa8M2/Y3nLg/v1kRDtUswJBLXhlAYULzgPw7q3LMnuQn4c1TrwCWXm477jdboPyDioTdb0P75rjie9YUm2OayMw6cplQR8O5eS3JsjpdvnUptN2Q95I4HQ0O9vkfOO0zJm9P8brNmM732bh6e42xJ9dpzvBAli0cacA798BHegIA8wG8KIS4AcAXiyGg7OJYae7ByQfGSonRttsYLyiIje1/arypSP3QuitnAIsvUbft8GbR/VBCNZOAUEV2N6w5OspjpCrKgU6u44Eu7dpCYS+E3Pocl9AxtjouX87+rJ6sU+05FZgqeoeP5mx7UFyrSHWmODbqZvMdyDNNr+jW5QpFsX5jZalwhBDPAXDGgkmmJxBCjAKw0hMAxqzmqPna/zw/Kwr/DWaKvTQcbgAWXAAsvChZ9nbLuUDrgoxtL+mwexzJu8t1If/9cKA39SPMNDs63hSODp2932+oGi8RC+Z6XFfL91eT/7AyudWvc5jHs2lPFznCDVXYqELi1ySeLWWpcDSo0hNMNV//HMB5RPQjAOrYLACI6DozDfW6w4dzSwHtNcWsX+T7wi0FdJL2E439ND6prww7dqGnsJ5+s8kCiUmzAADxivoMJ6aQY6MdzwyMFPY6+E217ayTFxny2lq2nm2p107PNi8f1+smV1cZpNdeTYG++5A+jC44bb4pSbfofKcnEEIMArjGw3m558NpWwT07s+4LuGVBIUQEGq3YZ1Sy3Zgsu93ieJgr7zonDrPWs/Jyhtv2geAtsWI7xlCatLpzkg2phQXiuUCWix0w1s8h8Udr2NmrmOrc3B2U2DZfBy5Pe/RDlLnzXaYFr2sSzVkGQPOJoMkahEsajavx0JSkgqn3NMTAMaO+A17jiLeeyCndlSZMfsihi7O+30oNeg036jCrQw4Zx4UwGiwCu/Vn+zSBwGRGgRInfsng1h5oapIm9yKxaBmBjiiMeV4cUrwOnHJ1Quy05E4r5AznCaXaBi2OrkKUeiQPCbZBjdN9VF8JpJJrUTSExhfY00khC+cNQvtfZt91d456YPY0nax+mCkFolTv4Qj1XOlnnLHunFtC5UuPxrLgSDtDCJsmPoZHK6Zn7HPTD+oKQ1ynvvsPqkuXIfftY3jEXnNzo1cL9n0Jru5qBTWcOQq1Y6gsEf63aNiOOtnjQdvuHL02CxLhWOmJ3gJwHwi2ktE1wghYgCs9ARbATw4XukJZB+sloHtAIDTZjVlrPX6lMtxqGYB+uUwME63zKg+M6GXqbwqbHly0b4iZT4IFXi/SyYzwUk5Zj8EgGXT1K7jOtPK6bPTN7cC/tNFlzpeLqen9UHkPuhNduT7yXUMdbrzZxc9QfIwrfDvKZcPrzJPa0V52KdUbMrStlDq6QmceL0BR0Kq7I5GXSs4pltLXswblQrPtpaaCDDvT4GqZuCFg6meMzRHlP0ekBFNFkILe9bG7BgcVa99tdapTSu6H/B8TcBQnfeQc9AbL7RfTZb3CWBPdAfkvr7gjBSeqwI7e749KG2uMxwnXm73fMQl8yJ2qybtuec+xkFflaXCKWkc3+IHZuojIVcECaPm7uBYwP3mSZqxXO4SLzufdx0ZSCub3VoDNLWDhAAgK5zs70hV8jU/A7Hcty5FroVuR7pucJCVmYxOYevcz3XeQ173e5QyOmV68gz7rDFXBeGsnesY6DUqtRtuHynh4QlrYR6cUrxc12mT1EnWShlWOIVAuinltZBoOGCLa7SssxG9Q2OojYawtu9EdVPmT9Awl6lcCFLosknKqHLPWOL6HTxyGRwy1ZXXX3YfGXQ9d2lHg7Jcdz2cIVgsVEoS0G+CXDxFNSPVD1hyGmMv8uRKIR5gnakmcp3hpMVf83APuu1LSW/Pv0xuMhzqzbyG49ws6gXnw6IXsWc0qR94dPdxKVCWazjlyvJpjbZBKkiExqoKhAIBXLC0XVnn7eaPYjDchIj5I1P9Fja3XYI9Dad6GricG9kAvSkp403v8sPsVAyuurS3KsLSonWmQU0nxlyNotANKB0+5AP0MyKdCVWXdjqiGUCdA4ff2GLadbgc3KWdNXOe4Tiqe3G1n+WyLyV9xpRftavK7eTEa7w0maZqp2kxcx3dOafMzLxePF6UriosW9zvlLpoGI1V4TTl0FQVwadP7URNJIRQIIDeJ2pwpH8EvdEpeKvlPKzCk0brirtsINKCgUiLUpk4UZl7nAu3yU+S4aYPkj5sQ7PCvNdW593zzB6zyl0OnYnM72Co24ugMwXqrrdu/UMXFVu3z8ip0KLhIMbiMa2Jz6nQopq9YL3D7ukg3HD2kes6QDYKy21/UaHDwHgxqTkl8FJnzKHIvChK3QNPOOTtGhQ6VJIKnuEUBPcbbMHkOsWCn0BrbRRVFSFUhAJoron4zm5ZV5n5+UGViEv3I830NFfhslmsQ2FWaJAWnDOF0pAHokwPlfmw27uhUyC6zb1OV18L3eA6vUltanPul2mvN67pB+epvemcM1VrY/BZc+3nZxWTzMT5JG59j25mLj8bbb3MBma5pC2ozyIywMJ2tWlUhZe5oXNmqdsnJeN8GPEit+6+PzqQfQr7QnNczHCIqBPAHTDir70thLh1nEXyhO0JxtMU236Syv1XZb7RPTH3DLrfuCNNi7C3P6TpJ73Npixzr2RyNvAaJj5bdLLq3Kt17uk6OXXlOhOpc49MTSSE/pFYWgZK6wnW+UARChBiWcUlSr9XrGuj8+QDvG+4BLzNBtyieDgfDrzcZl4e1CyioUBBQi3lc6d/S636vilGTLZMjL8EWZBFeoIlAB4SQlwNYHmBhcuunvOHdur1xh8MVwE/lgJrcdqZPx3w5ylmuXe2aUxusfkfw96GlbY9MxaqH7HsCOBHjkxBSgvthqybmejWSHSKqFojp07hOL2QEskcR/bzLMXkPN+KFuE0veQSekW3z8VNTfjxprJkc2YPlXFLM+586vdislOtN+qoqvD/jO4lzFWa00AO1q5JmnQGbmtfxaIsFQ78pyd4GcA1RPR7wFwMKST5yPhZ2Wj8mRC8b3qxTFcqP/1sbOaVFanbRPmUpGhS1Y/8FGe5ec/3YG6xYredOE0d8DOTQio2OhOl7to7E9BZe12cJlVL4aRlZjXfOm8Py5vcKU+nafLzY0pK9eX//vHz9G5NvNz2sniNggB4m+H4SZhoOXL4CeljrZGqzNkWzsSCutmpFxoce6UsmRdk8X3nm7JUOFmkJ7gKwDeFEB8G8LHiSeqR9qVA66LM52lY3tlgW8y1donXSF5OfjywLJt7e30UVRXBpNfLFR/oxOdPn2Hrd/GUOqzoTClG1SCmmoFYswAv6y/Wb7tcQnn4TczlPNu6Xs7BxlI4zvKUW7u9HWsNzjk4hgOyq31uxJMy6c9xDoBuJDQyy/iJ/OD2VVjHdN6DKizz4ErHfiQ3rNmLzmsS0Ju7slE4znQf1ufTeUPqLBiFYCKt4ajSE5xivn4SwLeI6FMAdusaIKLrAFwHAJ2dnYWRUsW884GA+ke0ckYjdu0btdks/mTJZGP28kYtMNKHVfNbsWp+a/L4is4GtNZGbPs+Llk+Nc0b6uozZqLfsRv/r86Zm3wdDQfxZ2fPTr53PqFFQkGcu9ge1PvcRW04Z2GrrezyD3TikBR5GjAWfs+Y05w2a/k/KzvSntQXT6nD+z1DOMWxiXbV/BZsfK8HTma31qT150ZDVVi5XjW5PqpMzzulIapMqjW7tQbvHOpPK18xvREb9qRHxp7XVottB/rSBsXFU+pw4NhwWjTweW212H1kMG0v1eKpdXi3ezBtrWTB5Fps3d9ri0sHpB5EdAOtavCz1omcWGObapZZVxlG79CYp4V2Z986Dzu/WMpfFVm9qSaCI30jvny1rNmCF49QC0t5qpRfZUVQuSZkXVedQ4mKJVPrsXnfsbTy9oZK7D82nPbgZ7nZ+1nDyhXKKsR8gfGSnoCIZgB4XAhxgvn+UgDnCyG+YL7/LIBThBBfykaGlStXinXr1vmv+NavgO53gNO/bLx/+v9lrnP21wE3M8FIP/Dij4zMnWd8JVUeGwHio0BkYoXbLxQjMeOH7TTxxOIJCKTPthIJYxUkbaaREEgIkbaGI4SAEN7TbccTAsNjcV+hUIQQOc/0hBB45/AAZrdUp7V1dGAUkXAgba1iaDSOobF4miNDPCGwfs9RrOhsSLsevcNjeL9nCAsmp896DxwbRiCQbvaNJwQ27e3B0o6GtOv+fs8Q9nQN4rTZ6ftMDhwbRu/wmNKzc9eRAbTURtIU7OBoDIf7RpRehW++34uW2kja2ooQAnu6BjG9qSrt2h0bHENciLRrNBKL48V3unDmnOa0e6xveAzdA6NKGY4OjKI2Gkq7rkOjccSFSPs8QgjEE+n3ZSIh0DUwqlwr3H6wD51NVXlzWiCi9UKIlbrjJTnDmQjpCYpCKGL8MZ7Q/ah0i/86xREIEAKK52Ii8rXYGwyQ77hb+TArEpE2eoIuaGdlRVBpygoGSBu+qS4aRt1k9UxAt54RDBCWd6rNVVMaKrW7+CfXR7Vt6kIQVVWEML1Jff1VDjeAce1maNpzumRbREJBfGh+q/JYbTSsXUNy+y50sqk2cQcCpHVMmetz60WulOUajoYSSU9QSEpvNsowDOOVslQ4JZ2eIJcQygzDMBOYkjSpZaLc0hPkTLjSUGSzPjTekjAMw2RNWSqc445AEFh1Y+bzGIZhSpiyNKmVNgTfay1lsr+EYRgmF1jhMAzDMEWBFQ7DMAxTFFjh5Bv2UmMYhlHCCodhGIYpChNS4RDRLCK6m4geMt9XE9HPiOgnRPTp8ZaPYRjmeKQsFI7f/DdmxOhrpFM/ASMfzrUALiqwtPDkpda6QKrCXmoMw0x8ykLhwH/+GycdSEWSzn+6Ph0JRVe1qpikDMMwE5+yUDhZ5L9xsheG0gFcPjMRXUdE64ho3eHDh3MVW820U1Kvg7zvlmGY44eyUDgaVPlvpgIAETUR0V0AlhPRTQB+DuCTRPRjAL/QNSiEWCOEWCmEWNnS0pKdVJm81GTz2Yorgdkfzq4fhmGYMqMkHrG95L/xgxCiC8D1juKrspEtJzK5R9e0GH8MwzDHASWhcDj/DcMwzMSnnE1qpZn/RggjE6eX8xiGYY4jykLhlHT+Gyf71rsftzJ0clpohmGOM0rCpJaJssx/o3KJDkWAxpnAwguBlgXpxxmGYSYwZaFwypJEHCDHBHL2hw0vtcknjI9MDMMw40jOJjUiKt5GynKnsnG8JWAYhhk38rGGw3FZtEiOAa0LgMbp4ycKwzDMOONZ4RDRXxDRKiJqdhxidysljsvSfuL4iMEwDFMi+Jnh3ANgEMDfE9HGgkgzkZk0a7wlYBiGGVf8KJx/AHA2gIcArCyMOPlBkZ7g42ZqggeI6NyiCXJ4W9G6YhiGKXX8KJzHAAwBOA3ADYURR02u6QmEEI+aqQmuB3BZUYQWAjj2XubzGIZhjhP8KJxFALYAuFsI8U8FkkfHPcgtPYHFN8w6hSNcWdDmGYZhypW8ruEQ0f+VXs/PSTKJXNMTkMFtAJ4QQmzQ9ZOX9ATTTzf+J2LZ1WcYhpmg+Nn4+V0A+2Gs4dhMakRUD+AHMELPDAHYBOAaFDZCsyo9wSmmPE0AbkEqPcEAgHMA1BPRHCHEXaoGhRBrAKwBgJUrV2bnfScSxv+tvwCi9Vk1wTAMMxHxo3COwFjDOdX8S5rVhBDHiOh3MOKdrQewBMDDXhsuUnqCO/y2kxWWwjm6G2hfWpQuGYZhygE/CuclAIsBhDXHDwH4GwCfhaGctgF43EvDEys9Ae+DZRiGUeFnDecSAJ0AXhZC3OI8KIT4DYBXhBBnAfgcgJr8iKilNNMTMAzDMEr8KJwRALcCcNvBWEdEK8xzq3MRTKas0hPIKaQ55w3DMEwSPya1XwH4DICZRPRNIcS3FefcAOCLMBTBk3mQD0CZpicAgJHe8ZaAYRimZMiocIjo5wC2A3gVwCNCCO1uRiHEGIq1OF8ODB0dbwkYhmFKBi8mtZ8D+E8AUQD/l4geJaJPF1ascoZNagzDMCq8mNQOAfg+DIVjeZ+dB+C+Aso1QWCFwzAMY5FxhmN6n70qeZ9VA/jzQgtWtrDTAMMwjBKvXmo27zMhRH8BZWIYhmEmIF4Vzg0AzgRwF9K9z0pup6MzPYFZVm3GSVtdPEl4hsMwDGPhSeEIIcaEEHcIIa4WQjzgOJaPNNWu5JqewOTrAB4stKzsNMAwDKOm4MoiT9yDHNITENFHAbwJwwGisPTK0XVY4TAMw1j42fg5bgghniOiGY7iZHoCACAiKz3Bm4omVsFwdlgEYIiIfiWEFWUzBRFdB+A6AOjs7MxO2ENbZcGza4NhGGYCUhYKR4Pn9ARCiJvN8isBHFEpGyBP6QkYhmEYJSWhcIqUngBCiHv8S5cLrLMYhmEsSkLhTKz0BBJsUmMYhklSLk4DKko/PYHacscwDHNcUhYKp6zSE8gk4uMtAcMwTMlQEia1TJRtegKGYRgmSVnMcMqWSW656hiGYY4vWOEUkqpJxv9geHzlYBiGKQFY4RSDQHC8JWAYhhl3WOEUEnaLZhiGScIKp6BYCqfkAmozDMMUnbLwUvMLEc0CcDOAeiHEpUQUAPBdAHUA1gkhflYUQawZDrHCYRiGKYsZTh7SE1wMIxLBGIyYa0WCZzgMwzAWZaFwkGN6AgDzAbwohLgBwBcLKKedaL3xv21x0bpkGIYpVcpC4QghngPQ7ShOpicQQowCsNITqNgL4Kj5unjb/yuqjf9TTypalwzDMKVKWSgcDar0BFMBIz0BEd0FMz0BgJ8DOI+IfgTgOV2DRHSdmYZ63eHDh3OXkNdwGIZhkpSE00CR0hM4U06r6nE+HIZhmAJREgpnwqYnYKcBhmGYJOVsUiv99AQMwzBMkrJQOGWbnuDIDuM/58VhGIYpDZNaJso2PUGXqXBiw+MrB8MwTAlQFjOcsof4MjMMw/BIWAxY4TAMw7DCKQqscBiGYVjhFAXe+MkwDMMKpyjwDIdhGGZiKhwimkVEdxPRQ+b7TiJ61Iw6fWOm+jlhBey0CcQZPxmGYcpC4eQhPcESAA8JIa4GsLygwkZq0st4hsMwDFMeCge5pyd4GcA1RPR7AE8WUE6gflp6GSschmGY8lA4eUhPcBWAbwohPgzgY4WTFGoHAVY4DMMw5aFwNPhJT/AkgK+YZbt1DeYnPQErHIZhGBUlEdqmSOkJLvVQL/f0BDzDYRiGUVISCmdCpSdo6ATwB3sZ78NhGIYpa5NaaaYniDakl7HCYRiGKQ+FU7bpCRiGYZgkJWFSy0TZpieI1gPDx8ZbCoZhmJKgLGY4ZUvTbCBcOd5SMAzDlASscApJ3wFev2EYhjEpC5Na2dL7/nhLwDAMUzLwDCfvZLd9h2EYZqLDCodhGIYpCqxwGIZhmKIwIddwiOjjMIJ01gG4G8bW/38DMArgGSHEfeMnHcMwzPFJWcxwssiH86gQ4loY8dQuA/AJGPlwrgVwUYGlLWzzDMMwZUpZKBxknw/nG+Y5HUhFlo4XVFKGYRhGSVkoHL/5cMjgNgBPCCE2wEhd0GHW037m/KQnYBiGYVSU8xqOKh/OKebrLwM4B0A9Ec0B8J8A7iSijwH4ha7BgqUnYBiGYUpD4RQgH84dAO5wFF+VjWwMwzBMfigJhTOh8uEwDMMwSspiDUdDaebDYS81hmEYJWWhcDgfDsMwTPlTEia1TJRtPhyGYRgmSVnMcMoK9lJjGIZRwgqHYRiGKQqscBiGYZiiUBZrOOUFm9QYpliMjY1h7969GB4eHm9Rjiui0Sg6OjoQDod91WOFwzBM2bJ3717U1tZixowZIF4/LQpCCHR1dWHv3r2YOXOmr7oTUuEo0hNUye+FEL8ZP+kYhskXw8PDrGyKDBGhqakJ2cSbLIs1nFzTEyjSFRQQTjHNMMWElU3xyfaal4XCQe7pCXTvGYZhmCJRFgon1/QEinQFSvKSniBYkV09hmHKlkcffRREhLfeekt5fNWqVVi3bl1RZfrHf/xH2/vTTz89+fprX/saFi9ejK997Ws4fPgwTjnlFCxfvhzPP/98QWUqC4WjQZWeYKr52kpPcCkRXa94r0QIsUYIsVIIsbKlpSU7qUKR7OoxDFO2rF27FmeeeSbWrl073qIkcSqcF198Mfl6zZo12LRpE26//Xb87ne/w5IlS/D666/jrLPOKqhMJeE0UKT0BM73DMNMIJ7ZdgiH+0by2mZLbQSr5re6ntPf348XXngBTz/9NC688EJ8+9vfxtDQEK666iq88cYbWLBgAYaGhpLnf/GLX8Rrr72GoaEhXHrppfj2t78NAJgxYwauuOIKPPHEEwiFQlizZg1uuukm7NixA1/72tdw/fXqZ+X9+/fjsssuQ29vL2KxGH784x/jl7/8JYaGhrBs2TIsXrwY9913H2pqatDf34+LLroI/f39OOmkk3DFFVfgX//1XzE0NIR169bhpZdeQmVlZf4uoIOSUDicnoBhmHLlf//3f3H++edj3rx5aGpqwvr16/Hss8+iqqoKW7duxaZNm7BixYrk+bfccgsmTZqEeDyOj3zkI9i0aROWLl0KAOjs7MTGjRvx1a9+FVdeeSX+8Ic/YHh4GCeccIJW4dx///0477zzcPPNNyMej2NwcBBnnXUW7rzzTmzcuDHt/Mceeww1NTXJY21tbVi3bh3uvPPOvF8bJyWhcLIkmZ4AhqK5HMCnxlckhmHGi0wzkUKxdu1a/OVf/iUA4PLLL8fatWuxY8cOfOUrXwEALF26NKlQAODBBx/EmjVrEIvFsH//frz55pvJ4xdddBEAYMmSJejv70dtbS1qa2sRiUTQ09ODhoaGtP5PPvlkXH311RgbG8PHP/5xLFu2rLAfOAfKQuGY6QlWAWgmor0AvimEuJuIrPQEQQA/5fQEDMMUk+7ubvz+97/H5s2bQUSIx+MgIixfvlx5/q5du/DP//zPeO2119DY2Igrr7zSFiUhEjHWgAOBQPK19T4Wiynb/OAHP4jnnnsOv/zlL3HllVfihhtuwOc+97k8fsr8URZOA0KIK4QQ7UKIsBCiQwhxt1n+KyHEPCHEbCHELeMtJ8MwxxcPPfQQPvvZz2LPnj3YvXs33nvvPcycORMnnXQS7r//fgDAli1bsGnTJgBAb28vqqurUV9fj4MHD+KJJ57IWYY9e/agra0N1157Lb7whS9gwwbDETccDmNsbCzn9vNJWcxwGIZhSpG1a9fi61//uq3sk5/8JF5//XUMDQ1h4cKFWLhwIU466SQAwIknnojly5djwYIFmDZtGs4444ycZXjmmWdw++23IxwOo6amBvfeey8A4LrrrsPSpUuxYsUK3HfffTn3kw9ICN4Zr2LlypUia7/5p/+f/f2HbspdIIZh0ti6dSsWLlw43mIcl6iuPRGtF0Ks1NUpC5MawzAMU/6wSY1hGKYM2Lx5Mz772c/ayiKRCF555ZVxksg/rHAYhmHKgCVLlij31ZQTE9KkRkQfJ6KfENEDRHSuWVZtxklbPd7yMQzDHI+UhcLJNT2BWfx1AA8WT2qGYRhGpiwUDnJMT0BEHwXwJoBDhReVYRiGUVEWCifX9AQwohScCiP0zbVEVBafm2EYZiJRzgOv5/QEQoibhRB/BeB+AD8RQiRUDeYlHw4ANM/Nvi7DMGVFTU2N6/F77rkHLS0tycjNl156KQYHBwEA3/rWtzB16lQsW7YMy5Ytw403GisDsVgMf/u3f4u5c+cmj91ySyqYyi233ILFixdj6dKlWLZsmW9PtUcffRRvvvlm8v3f//3f46mnngIAPP/881i8eDGWLVuGoaEhW+6cXCkJL7UipSeAEOKeDPXWAFgDGBs//fabZMml6Zs/GYY5brnsssuS0Zg/9alP4YEHHsBVV10FAPjqV7+Kv/7rv7ad/41vfAMHDhzA5s2bEY1G0dfXh3/5l38BALz00kt4/PHHsWHDBkQiERw5cgSjo6O+5Hn00UexevVqLFpkrEJ85zvfSR677777cNNNN+Ezn/kMACN3Tnd3N4LBYHYfXqIkFA6nJ2AYJme2PwX0H8xvmzVtwNxshic1sVgMAwMDaGxs1J4zODiIn/zkJ9i9ezei0SgAoLa2Ft/61rcAGPlvmpubk8E9m5ubXfu88cYb8dhjjyEUCuHcc8/FJz7xCTz22GN49tln8Q//8A94+OGH8d3vfherV69GT08PHnzwQfz617/GE088gb6+vmTunJtuugmXXXaZa1+ZKAmFkyWcnoBhmLLggQcewAsvvID9+/dj3rx5uPDCC5PHvv/97+O//uu/AAC33XYb2tvb0dnZidraWmVb5557Lr7zne9g3rx5OOecc3DZZZfh7LPPVp7b1dWFRx55BG+99RaIKJni4KKLLsLq1atx6aWX2s7/whe+gBdeeMF2TM6dkytloXDKMj1B22LgYOmIwzATnjzORPKNZVITQuAv/uIvcPvttyfXa5wmNSuytMV//Md/4Ic//CG6urrw4osvYtq0aVi/fj2ef/55PP3007jssstw66234sorr0zrt76+HtFoFNdccw1Wr16N1avHdxtiWTgNlGV6gnnnZz6HYZgJx3vvvZdc6L/rrrtsx4gIF154IZ577jlt/Tlz5uDdd99FX18fAOCqq67Cxo0bUV9fj3g8DgAIBoNYtWoVvv3tb+POO+/Eww8/rGwrFArh1VdfxaWXXorHH38c558/vuNSWcxwypJA7gtsDMOUH9OmTbOZoO655x7b8RdeeAGzZ8/W1q+qqsI111yDL33pS/j3f/93RKNRxOPxpGPAtm3bEAgEMHeu4Q27ceNGTJ8+XdlWf38/BgcHccEFF+CMM87ArFmzABhrQpZCKyascAoNKx6GmdAMDg6io6Mj+f6GG27ADTfcYDvHWsNJJBLo6OhIU0JObrnlFvzd3/0dTjjhBNTW1qKyshKf//znMWXKFGzevBlf/vKX0dPTg1AohDlz5mDNmjXKdvr6+nDxxRdjeHgYQgh873vfA2Ckwr722mtxxx134KGHHsrtAviA8+FoyCkfDgAk4sCz/wSEIsBZN2Q+n2EY33A+nPEjm3w4PMMpFIEgMPvDQJN+6swwDHM8wQqnkHSeMt4SMAxznHDJJZdg165dtrLbbrsN55133jhJlM6EVDhE9HEAHwNQB+BuAE8B+K75fp0Q4mfjJx3DMEz+eeSRR8ZbhIyUhVt0HtITXAwjEsEYjJhrDMNMEHgduvhke83LQuEgx/QEAOYDeFEIcQOALxZcWoZhikI0GkVXVxcrnSIihEBXV1cy7I4fysKkJoR4johmOIqT6QkAgIis9ARvEhEBuBVmegJTEVnR7eJFEpthmALT0dGBvXv3Iqfo7oxvotGozRXcK2WhcDSo0hNYq/RWeoJ6IpoD4F4APyKiswBot/gS0XUArgOAzs7OQsjMMEweCYfDmDlz5niLwXikJBROkdITXOOhXn7SEzAMwzBplITC4fQEDMMwE59ycRpQkUxPQEQVMNITPDbOMjEMwzAayiK0jZyeAMBBpNITXADgB0ilJ8hbxGgiOgxgT5bVmwEcyZcseYTl8gfL5Q+Wyx8TUa7pQogW3cGyUDjlBhGtc4snNF6wXP5gufzBcvnjeJSrnE1qDMMwTBnBCodhGIYpCqxwCoM6OcX4w3L5g+XyB8vlj+NOLl7DYRiGYYoCz3AYhmGYosAKJ4/oolcXsL9pRPQ0Eb1JRH8kor80y79FRPuIaKP5d4FU5yZTvm1EdJ5UnlfZiWg3EW02+19nlk0iot8S0Xbzf6NZTkR0h9n3JiJaIbXzefP87UT0+Rxlmi9dk41E1EtEfzUe10sVAT2f14eITjKv/w6zLuUg1+1E9JbZ9yNE1GCWzyCiIem63ZWpf91nzFKuvH1vZOzne8Usf4CMvX3ZyvWAJNNuIto4DtdLNzaM7z0mhOC/PPzB2Av0DoBZACoAvAFgUYH7bAewwnxdC+BtGJGzvwXgrxXnLzLligCYacobLITsAHYDaHaU/ROAG83XNwK4zXx9AYAnABCAUwG8YpZPArDT/N9ovm7M4/d1AMD08bheAD4IYAWALYW4PgBeNc8ls+6f5CDXuQBC5uvbJLlmyOc52lH2r/uMWcqVt+8NwIMALjdf3wXgi9nK5Tj+LwD+fhyul25sGNd7jGc4+SMZvVoIMQrAil5dMIQQ+4UQG8zXfQC2wghqquNiAP8thBgRQuwCsMOUu1iyXwzASn73MwAfl8rvFQYvA2ggonYA5wH4rRCiWwhxFMBv4UhTkQMfAfCOEMJtc2/BrpcQ4jkA3Yr+cr4+5rE6IcTLwhgZ7pXa8i2XEOI3QoiY+fZlGGGktGToX/cZfcvlgq/vzXwy/zCAh/Ipl9nunwJY69ZGga6XbmwY13uMFU7+UEWvdhv88woZ6RuWA3jFLPqSOTX+qTQN18lYCNkFgN8Q0XoyonADQJsQYr/5+gCAtnGQy+Jy2AeC8b5eQP6uz1TYEw3m87pdDeNp1mImEb1ORM+SEY3dklfXv+4zZks+vrcmAD2SUs3X9ToLwEEhxHaprOjXyzE2jOs9xgpnAkBENQAeBvBXQoheAD8GMBvAMgD7YUzri82ZQogVMBLk/QURfVA+aD4VjYuLpGmfvwjA/5hFpXC9bIzn9dFBRDcDiAG4zyzaD6BTCLEcwA0A7ieiOq/t5eEzltz35uAK2B9qin69FGNDTu3lCiuc/DEu0auJKAzjhrpPCPFzABBCHBRCxIUQCQA/gWFKcJMx77ILIfaZ/w8BeMSU4aA5FbfMCIeKLZfJnwDYIIQ4aMo47tfLJF/XZx/sZq+c5SOiKwGsBvBpc6CCabLqMl+vh7E+Mi9D/7rP6Js8fm9dMExIIUd51phtfQLAA5K8Rb1eqrHBpb2i3GOscPJH0aNXmzbiuwFsFUJ8Typvl067BIDlQfMYgMuJKEJEMwHMhbHwl1fZiaiaiGqt1zAWnbeYbVpeLp8HYOU6egzA50xPmVMBHDOn/b8GcC4RNZrmknPNslyxPXmO9/WSyMv1MY/1EtGp5j3yOakt3xDR+QD+BsBFQohBqbyFjFTvIKJZMK7Pzgz96z5jNnLl5XszFejTAC7Nh1wm5wB4SwiRNDsV83rpxgaX9opzj2XyKuA/738wPD3ehvHkcnMR+jsTxpR4E4CN5t8FAP4TwGaz/DEA7VKdm035tkHyKsmn7DC8gN4w//5otQfDVv47ANsBPAVgkllOAP7V7HszgJVSW1fDWPTdAeCqPFyzahhPtPVSWdGvFwyFtx/AGAz79zX5vD4AVsIYgN8BcCfMTd5ZyrUDhh3fusfuMs/9pPn9bgSwAcCFmfrXfcYs5crb92bes6+an/V/AESylcssvwfA9Y5zi3m9dGPDuN5jHGmAYRiGKQpsUmMYhmGKAischmEYpiiwwmEYhmGKAischmEYpiiwwmEYhmGKAischikARPSi+X8GEX0qz23/raovhil12C2aYQoIEa2CEdF4tY86IZGK66U63i+EqMmDeAxTVHiGwzAFgIj6zZe3AjiLjPwnXyWiIBn5ZV4zg07+mXn+KiJ6nogeA/CmWfYoGcFP/0hmAFQiuhVApdnefXJf5i7x24loCxl5Si6T2n6GiB4iI6/NfebucBDRrWTkTNlERP9czGvEHH+EMp/CMEwO3AhphmMqjmNCiJOJKALgD0T0G/PcFQBOEEZIfQC4WgjRTUSVAF4jooeFEDcS0ZeEEMsUfX0CRiDLEwE0m3WeM48tB7AYwPsA/gDgDCLaCiMkzAIhhCAzsRrDFAqe4TBMcTkXRsyqjTDCxTfBiKkFAK9KygYAvkJEb8DIQTNNOk/HmQDWCiOg5UEAzwI4WWp7rzACXW6EkQzsGIBhAHcT0ScADKY3yTD5gxUOwxQXAvBlIcQy82+mEMKa4QwkTzLWfs4BcJoQ4kQArwOI5tDviPQ6DiODZwxGhOWHYESCfjKH9hkmI6xwGKaw9MFI8WvxawBfJCN0PIhonhlR20k9gKNCiEEiWgAjla/FmFXfwfMALjPXiVpgpD9+VScYGblS6oUQvwLwVRimOIYpGLyGwzCFZROAuGkauwfAD2GYszaYC/eHoU7N+ySA6811lm0wzGoWawBsIqINQohPS+WPADgNRpRuAeBvhBAHTIWlohbA/xJRFMbM64asPiHDeITdohmGYZiiwCY1hmEYpiiwwmEYhmGKAischmEYpiiwwmEYhmGKAischmEYpiiwwmEYhmGKAischmEYpiiwwmEYhmGKwv8HP9mRpikARm0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plt.plot(s1, alpha = 0.8, label = 'Adam_non_stiff')\n",
    "plt.plot(s2, alpha = 0.5, label = 'Adam_stiff')\n",
    "# plt.plot(s3[:,0], alpha = 0.5, label = 'L-BFGS_non_stiff')\n",
    "plt.plot(s4[:,0], alpha = 0.5, label = 'L-BFGS_stiff')\n",
    "# plt.hlines(np.mean(s1), 0, 4999, colors='blue', linestyles='dashed')\n",
    "# plt.hlines(np.mean(s2), 0, 4999, colors='orange', linestyles='dashed')\n",
    "plt.yscale('symlog')\n",
    "\n",
    "plt.xlabel('iterations')\n",
    "plt.ylabel(r'$\\frac{d J_{PINN}}{dt}$')\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig('prove_stiffness/compare_DJ_dt.png', dpi= 500)\n",
    "\n",
    "print(np.mean(s1))\n",
    "print(np.mean(s2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
